The first integral can be written as
$$\int dU \, \operatorname{tr}(U \rho U^\dagger M^\dagger M)
= \operatorname{tr}\left[\left(
\int dU\, U\rho U^\dagger \right)M^\dagger M\right],$$
where the integration is performed over the uniform Haar measure in the space of unitaries. You then get the result using the identity
$$
\int dU\, UX U^\dagger = \frac{\operatorname{tr}(X) I}{d}$$
for any linear operator $X$, for the special case of $\rho$ pure. This has been discussed and proved multiple times on the site already, see e.g. Density matrices of multiples copies of a single Haar-Random state and links therein.
For the other identity, write the LHS as
$$\int dU \operatorname{tr}(U\rho U^\dagger M)^2
= \sum\int dU\, U_{ij} U_{k\ell} \bar U_{mn} \bar U_{pq} \rho_{jn} \rho_{\ell q} M_{mi} M_{pk}$$
where I wrote expliciting the matrix components. Using the formulae to integrate polynomials of the unitary matrices, which you can find e.g. in the Wikipedia page, you have
$$
\int_{U_d} dU U_{ij} U_{k\ell} \bar U_{mn}
\bar U_{pq}= \frac{1}{d^2-1}\left[ (\delta_{im}\delta_{jn} \delta_{kp}\delta_{\ell q} + \delta_{ip}\delta_{jq} \delta_{km}\delta_{\ell n} ) - \frac1d
(\delta_{im} \delta_{jq} \delta_{kp}\delta_{\ell n}
+\delta_{ip} \delta_{jn} \delta_{km}\delta_{\ell q})
\right]$$
which contracting the indices should give back your expression.
