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 can be either reduced to Goncharov’s polylogarithms or 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.