We calculate three-loop photon spectral density in QED with N different species of electrons. The obtained results were expressed in terms of iterated integrals, which are either  reduce to Goncharov's polylogarithms or can be written in terms of one-fold integrals of harmonic polylogarithms and complete elliptic integrals. In addition we provide threshold and high-energy asymptotics of the calculated spectral density. It is shown, that the use of the obtained spectral density correctly reproduces separately calculated moments of corresponding photon polarization operator.