The general formula in the circular ensembles follows from equation 2.10 of arXiv:cond-mat/9612179:
$$P(H)\propto\text{det}\,(1-H)^{\tfrac{1}{2}\beta(N-2M+1-2/\beta)},\;\;N-2M\geq 0,\qquad\qquad(\ast)$$
where $H=AA^\dagger$ and $A$ is an $M\times M$ principal submatrix of an $N\times N$ unitary matrix $U$ (circular unitary ensemble, $\beta=2$) or unitary symmetric matrix (circular orthogonal ensemble, $\beta=1$) or unitary selfdual matrix (circular symplectic ensemble, $\beta=4$).
The case in the OP is $\beta=1$, in that case $H=AA^\dagger=AA^\ast$ (in the physics notation that $\ast$ is complex conjugation and $\dagger$ is Hermitian conjugate). The assumption that $N-2M\geq 0$ is needed, because if $M>N/2$ there will be an additional set of $M-N/2$ eigenvalues of $H$ that are pinned to unity.
The case that $U$ is real orthogonal (circular real ensemble, $H=AA^\dagger=AA^\top$) has an additional factor,
$$P(H)\propto\text{det}\,H^{-1/2}\,\text{det}\,(1-H)^{\tfrac{1}{2}(N-2M-1)}.\qquad\qquad(\ast\ast)$$ That case is considered in an earlier MO posting. Note that $(\ast\ast)$ agrees with the corresponding result in the OP (the second equation), which gives $P(A)$ instead of $P(H)$.
To make contact with the physics literature: $U$ represents a scattering matrix, $A$ is the reflection matrix, and the eigenvalues of $1-H=1-AA^\dagger$ are the transmission eigenvalues $T_1,T_2,\ldots T_M$. Equation (2.10) in the cited paper gives the joint probability distribution of the $T_i$'s, which can equivalently be expressed as the equation $(\ast)$ above.
