We calculate the non-singlet, the pure singlet contribution, and their interference term, at O(α2) due to electron-pair initial state radiation to e+e− annihilation into a neutral vector boson in a direct analytic computation without any approximation. The correction is represented in terms of iterated incomplete elliptic integrals. Performing the limit s≫me2 we find discrepancies with the earlier results of Ref. [1] and confirm results obtained in Ref. [2] where the effective method of massive operator matrix elements has been used, which works for all but the power corrections in m2/s. In this way, we also confirm the validity of the factorization of massive partons in the Drell–Yan process. We also add non-logarithmic terms at O(α2) which have not been considered in [1]. The corrections are of central importance for precision analyses in e+e− annihilation into γ⁎/Z⁎ at high luminosity. Keywords: Higher order QED correction, e+e− annihilation, Precision physics at ILC, FCC, Iterated integrals