Integral around unit sphere of inner product

The integral is linear in $M$, so we only have to calculate it for canonical matrices $A_{kl}$ spanning the space of matrices, with $(A_{kl})_{ij}=\delta_{ik}\delta_{jl}$. The integral vanishes by symmetry for $k\neq l$, since for every point on the sphere with coordinates $x_k$ and $x_l$ there's one with $x_k$ and $-x_l$.

So we only have to calculate the integral for $k=l$. By symmetry, this is independent of $k$, so it's just $1/n$ of the integral for $M$ the identity. But that's just the integral over $1$, which is the surface area $A$ of the sphere.

Then by linearity the integral for arbitrary $M$ is the sum of the diagonal elements, i.e. the trace, times the coefficient $A/n$.

As a function of $M$ your integral is linear, and is invariant under conjugation by orthogonal transformations ($C_R: M \mapsto R^{T} M R$). Now the average of $C_R(M)$ over all orthogonal transformations $R$ (using Haar measure) is $\text{Tr}(M) I/n$ (it must be invariant under all $C_R$ so it is a multiple of $I$, and the trace is preserved). So the integral is the same as it would be for $\text{Tr}(M)I/n$, which is $\text{Tr}(M)/n$ times the area of the sphere.