  We introduce novel relations between the derivatives
[$\partial^{n}\rho(\lambda, m_{l})/\partial m_{l}^n$] of the Dirac eigenvalue
spectrum [$\rho(\lambda, m_{l})$] with respect to the light sea quark mass
($m_{l}$) and the $(n+1)$-point correlations among the eigenvalues ($\lambda$)
of the massless Dirac operator. Using these relations we present lattice QCD
results for $\partial^{n}\rho(\lambda, m_{l})/\partial m_{l}^n$ ($n=1, 2, 3$)
for $m_{l}$ corresponding to pion masses $m_\pi=160-55$ MeV, and at a
temperature of about 1.6 times the chiral phase transition temperature.
Calculations were carried out using (2+1) flavors of highly improved staggered
quarks with the physical value of strange quark mass, three lattice spacings
$a=0.12, 0.08, 0.06$ fm, and lattices having aspect ratios $4-9$. We find that
$\rho(\lambda\to0, m_{l})$ develops a peaked structure. This peaked structure
arises due to non-Poisson correlations within the infrared part of the Dirac
eigenvalue spectrum, becomes sharper as $a\to0$, and its amplitude is
proportional to $m_{l}^2$. We demonstrate that this $\rho(\lambda\to0,m_l)$ is
responsible for the manifestations of axial anomaly in two-point correlation
functions of light scalar and pseudoscalar mesons. After continuum and chiral
extrapolations we find that axial anomaly remains manifested in two-point
correlation functions of scalar and pseudoscalar mesons in the chiral limit.
