We investigate general properties of the eigenvalue spectrum for improved staggered quarks. We introduce a new chirality operator $[gamma_5 otimes 1]$ and a new shift operator $[1 otimes xi_5]$, which respect the same recursion relation as the $gamma_5$ operator in the continuum. Then we show that matrix elements of the chirality operator sandwiched between two eigenstates of the staggered Dirac operator are related to those of the shift operator by the Ward identity of the conserved $U(1)_A$ symmetry of staggered fermion actions. We perform a numerical study in quenched QCD using HYP staggered quarks to demonstrate the Ward identity. We introduce a new concept of leakage patterns which collectively represent the matrix elements of the chirality operator and the shift operator sandwiched between two eigenstates of the staggered Dirac operator. The leakage pattern provides a new method to identify zero modes and non-zero modes in the Dirac eigenvalue spectrum. This method is as robust as the spectral flow method but requires much less computing power. Analysis using a machine learning technique confirms that the leakage pattern is universal, since the staggered Dirac eigenmodes on normal gauge configurations respect it. In addition, the leakage pattern can be used to determine a ratio of renormalization factors as a by-product. We conclude that it might be possible and realistic to measure the topological charge $Q$ using the Atiya-Singer index theorem and the leakage pattern of the chirality operator in the staggered fermion formalism.