The following argument is quite similar to Carlo Beenakker's.
For simplicity, I only consider the real case. I'm also going to use different symbols for the sizes of the matrices. Let $X$ be an $n \times d$ random matrix with iid rows from $N(0,\Sigma)$ let $Y$ be a $m \times d$ random matrix independent of $X$, with iid rows from $N(0,S)$. Set $\widehat \Sigma := X^\top X/n$ and $\widehat S := Y^\top Y/m$. Let $F(a)$ be the expression you wish to evaluate and set $\lambda=1/(na)$. Note that
$$
F(a) = \frac{m}{na} \mathbb E\,[\operatorname{tr} \widehat S(\widehat\Sigma + \lambda I_d)^{-1} ]
$$
Then, in the limit $n,d \to \infty$ such that $d/n \to \phi \in (0,\infty)$, we have
$$
\mathbb E\,[\operatorname{tr} \widehat S (\widehat\Sigma + \lambda I_d)^{-1} \mid Y] \simeq \operatorname{tr} \widehat S(\Sigma + \kappa I_d)^{-1},
$$
where $\kappa = \kappa(\lambda,n)$ is the unique nonnegative solution to the equation $\kappa - \lambda = \kappa \operatorname{df}_1(\kappa)/n$, where
$$
\operatorname{df}_1(\kappa) := \operatorname{tr}\Sigma (\Sigma + \kappa I_d)^{-1}.
$$
We deduce that
$$
\mathbb E\,[\operatorname{tr} \widehat S(\widehat \Sigma + \lambda I_d)^{-1} ] \simeq \mathbb E\,[\operatorname{tr} \widehat S (\Sigma + \kappa I_d)^{-1}] = \mathbb E\,[\operatorname{tr}S (\Sigma + \kappa I_d)^{-1}].
$$
In your question, you have $\Sigma=S=I_d$, and so
$$
\mathbb E\,[\operatorname{tr}\widehat S (\widehat \Sigma + \lambda I_d)^{-1}] \simeq \frac{d}{1+\kappa}.
$$
Moreover, $\kappa$ is now given by the following well-known formula which is reminiscent of the Marchenko-Pastur distribution
$$
\kappa = \frac{\lambda + \overline \phi + \sqrt{(\lambda + \overline\phi)^2 + 4\lambda}}{2},
$$
where $\overline \phi := \phi-1$.
Simplify...