Let $X=A+\mathrm{i}B$, with $A,B\in\mathbb{R}^{N\times M}$ whose entries are standard normal centered gaussian variables.
We thus have:
$$X^\dagger X=(A^\top-\mathrm{i}B^\top)(A+\mathrm{i}B)=A^\top A+B^\top B+\mathrm{i}(A^\top B-B^\top A)$$
We're interested in:
$$\begin{align}\mathbb{E}\left[\frac{X^{\dagger}X}{\mathrm{tr}\left(X^\dagger X\right)}\right]&=\mathbb{E}\left[\frac{A^\top A+B^\top B+\mathrm{i}(A^\top B-B^\top A)}{\mathrm{tr}\left(A^\top A\right)+\mathrm{tr}\left(B^\top B\right)+\mathrm{i}\underbrace{\left(\mathrm{tr}\left(A^\top B\right)-\mathrm{tr}\left(B^\top A\right)\right)}_{0}}\right]
\end{align}$$
Let us perform an analysis on a per-component basis. On the diagonal, say for the $x^{\text{th}}$ line and column, the coefficient in this density matrix is the expected value of:
$$\frac{1}{\sum\limits_{j=1}^M\sum\limits_{i=1}^N\left(A_{i,j}^2+B_{i,j}^2\right)}\left[\sum_{i=1}^{N}A_{i,x}^2+\sum_{i=1}^{N}B_{i,x}^2+\mathrm{i}\underbrace{\left(\sum_{i=1}^NA_{i,x}B_{i, x}-\sum_{i=1}^NB_{i,x}A_{i, x}\right)}_0\right]$$
where all the $A_{i,j}$ and $B_{i,j}$ are centered standard normal random variables. We're thus interested in:
$$\mathbb{E}\left[\frac{\sum\limits_{i=1}^NA_{i,x}^2+\sum\limits_{i=1}^NB_{i,x}^2}{\sum\limits_{j=1}^M\sum\limits_{i=1}^NA_{i,j}^2+B_{i,j}^2}\right]$$
We can write it as:
$$\frac{S_1}{S_1+S_2}$$
where $S_1\sim\chi^2(2N)$ and $S_2\sim\chi^2(2MN-2N)$. Because of the link between the $\chi^2$ and $\Gamma$ distribution, we have $S_1\sim\Gamma(N, 2)$ and $S_2\sim\Gamma(N(M-1), 2)$, which means that $\frac{S_1}{S_1+S_2}\sim\beta(N, N(M-1))$, whose expected value is $\frac{N}{N(M-1)+N}=\frac{1}{M}$. Thus, every diagonal coefficient in this density matrix is equal to $\frac1M$, which was expected: there's no reason to prefer a computational state to another.
If we now take the coefficient at line $x$ and column $y$, the expression of this coefficient is:
$$\frac{1}{\sum\limits_{j=1}^M\sum\limits_{i=1}^N\left(A_{i,j}^2+B_{i,j}^2\right)}\left[\sum_{i=1}^{N}A_{i,x}A_{i,y}+\sum_{i=1}^{N}B_{i,x}B_{i,y}+\mathrm{i}\left(\sum_{i=1}^NA_{i,x}B_{i, y}-\sum_{i=1}^NB_{i,x}A_{i, y}\right)\right]$$
We thus want to compute:
$$\mathbb{E}\left[\frac{A_{i,x}A_{i,y}}{\sum\limits_{j=1}^M\sum\limits_{i=1}^NA_{i,j}^2+B_{i,j}^2}\right]=\mathbb{E}\left[\frac{XY}{X^2+Y^2+S}\right]$$
since this will give us the expectation we're looking for by linearity. Here, $X$ and $Y$ are independent and sampled from a centered standard gaussian, while $S\sim\chi^2(NM-2)$ is also independent from $X$ and $Y$. What is important here is that $X$ (and $Y$, though it doesn't matter because of the independence) is centered around $0$. Indeed, computing this expectation yields, using the fact that $f_{X,Y,S}(x,y,s)=f_X(x)f_Y(y)f_S(s)$:
$$\mathbb{E}\left[\frac{XY}{X^2+Y^2+S}\right]=\int_{0}^{+\infty}\left[\int_{-\infty}^{+\infty}\left(y\int_{-\infty}^{+\infty}\frac{x}{x^2+y^2+s}f_X(x)\,\mathrm{d}x\right)f_Y(y)\,\mathrm{d}y\right]f_S(s)\,\mathrm{d}s$$
Since $f_X$ is even, clearly the middle integrand is odd, which means that the integral is nil. All in all, this expectaion is equal to $0$.
Thus, the final density matrix is the maximally mixed state $\frac1MI_M$.
Now that we know that $\rho$ is the maximally mixed state, it is easy to compute $F(\rho, \sigma)$:
$$F(\rho, \sigma)=\frac1M\mathrm{tr}(\sqrt{\sigma})=\frac1M\sum_{i=1}^M\sqrt{\lambda_i}$$
where the $\lambda_i$ are the eigenvalues of $\sigma$.
For the average Fidelity between $\sigma$ and some $\rho$ taken from $\mu$, I'm not sure much can be said (maybe someone will add to this answer). If $\sigma=|\psi\rangle\langle\psi|$ is a pure state, this will give $\int\langle\psi|\rho|\psi\rangle\,\mathrm{d}\mu(\rho)=\langle\psi|\int\rho\,\mathrm{d}\mu(\rho)|\psi\rangle=\frac1M$, which corresponds to the previous case, but I'm not sure this is still true for general mixed states.
