The formula
$\mathbb E[(\sum_{i=1}^n x_ix_i^T)^{-1}] = I_p/(n-p-1)$ is exact. It is known on Wikipedia as the expectation of the inverse wishart distribtion: https://en.wikipedia.org/wiki/Inverse-Wishart_distribution
To prove it, outside of the diagonal for some $i\ne j$: with $X$ the matrix with rows $x_1,...,x_n$ and $e_j$ the $j$-th canonical vector, the fact that $X(I_p - 2e_je_j^T)$ has the same distribution as $X$ (i.e., changing the sign of the $j$-th column) gives that the off-diagional $i,j$-th element is 0.
For the diagonal elements, by symmetry it is enough to compute $trace[(X^TX)^{-1}]$ as asked in the question.
By Stein's formula $\mathbb E[x_i^T F(X)]=\sum_j \mathbb E[ e_j^T \frac{\partial}{\partial x_{ij}} F(X)]$ for $F$ valued in $R^p$ we find
\begin{align}
\mathbb E[x_i^T(X^TX)^{-1}x_i]
=\mathbb E\Big[\sum_{j=1}^p e_j^T(X^TX)^{-1} e_j - e_j^T(X^TX)^{-1}(e_j x_i^T + x_i e_j^T)(X^TX)^{-1} x_i\Big].
\end{align}
Summing over $i=1,...,n$
$$\mathbb E\sum_{i=1}^n x_i^T(X^TX)^{-1}x_i
= n \mathbb E\Big[Tr[(X^TX)^{-1}]
- Tr[(X^TX)^{-1}]\sum_{i=1}^n x_i^T(X^TX)^{-1}x_i
- Tr[(X^TX)^{-1}\sum_{i=1}^nx_ix_i^T(X^TX)^{-1}]
\Big]
$$
or equivalently $p=(n-p-1)\mathbb E Tr[(X^TX)^{-1}]$ which is the desired identity.