This is easy to do using a graphical calculus for contractions of old-fashioned tensors.
Here is a proof. When I have time I will try to make nicer pictures.
Represent a matrix element $M_{ij}$ by a triangular box labelled $M$ and two strands coming out carrying the indices $i$ and $j$, a bit like $i-<M|-j$ (sorry can't complete the two missing corners of the triangle), then see that picture rotated clockwise by 90 degrees.
Then the determinant is given by $$ {\rm det}(M)= \ \ \begin{array}{ccccccc} & & --& -- &-- & & \\ & / & --&-- & --&\backslash & \\ / & / & & & & \backslash & \backslash \\ | & | & & & & | & | \\ | & | & & & & M & M \\ | & | & & & & | & | \\ | & | & & & & = & = \\ | & | & & & & | & | \\ \backslash & \backslash & & & &/ & / \\ & \backslash & --& -- &-- &/ & \\ & &-- & -- &-- & & \end{array} $$
This picture is for $n=2$. The row of equal signs represents a single box for a normalized antisymmetrizer with $n$ strands going through it. There is also a "Markov trace" or closing the loops without crossings.
Now take $M=X^2+Y^2$ and expand by multilinearity. That means that, in the previous picture, each $$ \begin{array}{c} | \\ M \\ | \end{array} $$ becomes $$ \begin{array}{c} | \\ X \\ | \\ X \\ | \end{array} $$ or $$ \begin{array}{c} | \\ Y \\ | \\ Y \\ | \end{array} $$ Now the placement of the pairs of $X$'s and $Y$'s do not matter, because exchanging positions creates two twists, one on top and one below the antisymmetrizer. Undoing them gives a factor $(-1)^2$, i.e., does nothing. So we get a sum $$ \sum_{n_1+n_2=n}\frac{n!}{n_1! n_2!}\times $$ the expectation of the picture above with the $M$'s (now $n$ of them) are replaced by something like $$ \begin{array}{} | & | & \cdots &|&|&|& \cdots&|& \\ X& X &\cdots& X &Y& Y&\cdots& Y \\ | & | & \cdots &|&|& |&\cdots&|& \\ X& X &\cdots& X &Y& Y&\cdots& Y \\ | & | & \cdots &|&|&|&\cdots&|& \end{array} $$ where there $n_1$ vertical strands of $X$'s followed by $n_2$ strands of $Y$'s.
Now we apply the Isserlis-Wick Theorem to perform the Gaussian integrals as a sum over perfect matchings of $X$'s among themselves and of $Y$'s among themselves. The main ingredient is the graphical identity $$ \mathbb{E}\left[ \begin{array}{c} | \\ X \\ | \end{array} \begin{array}{c} | \\ X \\ | \end{array} \right]= \begin{array}{} | & & | \\ - & -- & - \\ & & \\ - & -- & - \\ | & & | \end{array} $$ namely, with indices, $$ \mathbb{E}[X_{ij}X_{k\ell}]=\delta_{ik}\delta_{j\ell} $$ A better picture than using "$--$"'s and "$|$" would show a cup $\bigcup$ on top and a cap $\bigcap$ on the bottom.
Notice that two $X$'s in the same row cannot contract, because otherwise we would get a symmetric Kronecker delta contracted directly with an antisymmetrizer and this gives zero. So the $i$-th $X$ in the bottom row should contract with the $\sigma_1(i)$-th $X$ in the top row and the $j$-th $Y$ in the bottom row should contract with the $\sigma_2(j)$-th $Y$ in the top row. Summing over Wick contractions thus reduces to summing over two permutations $\sigma_1\in S_{n_1}$ and $\sigma_2\in S_{n_2}$.
Now the hard part, and sorry the pictures get a bit tricky (try to draw a cycle of length four), the contribution of the permutations can be computed as follows. After straightening or pulling on the strands, one gets the compositional squares of the two permutations. In the Schur-Weyl picture we would have something like $$ \omega= (\sigma_1\circ \sigma_1) \otimes (\sigma_2\circ\sigma_2) $$ as an endomorphism of $\wedge^n(\mathbb{R}^n)$.
Another way to see this is as follows. Consider the $i$-th "$|$" below a bottom row $X$. It will follow a cap going to the strand "$|$" right below the $\sigma_1(i)$-th top row $X$. You have to go down that strand and then follow the cup that takes you to the "$|$" sitting on top of the $\sigma_1(\sigma_1(i))$-th top row $X$. At the end of the day after pulling on the strands, eliminating curlicues $X|$ and zigzags $/|/$, you have that the $i$-th strand at the very bottom connects to the $\sigma_1(\sigma_1(i))$-th strand at the very top. That was for the $X$ compartment, but the same is done independently in the $Y$ compartment.
Then the evaluation of the the contribution of $\omega$ is the sign of the permutation $\omega$, i.e., $1$ times the trace of the identity on $\wedge^n(\mathbb{R}^n)$ which is also $1$.
In sum, and for the general case of $m$ matrices, we have $$ \mathbb{E}[{\rm det}(X_1^2+\cdots+X_m^2)]= \sum_{n_1+\cdots+n_m=n} \binom{n}{n_1,\ldots,n_m} n_1!\cdots n_m!\ \times 1 $$ $$ =n!\sum_{n_1+\cdots+n_m=n} 1=n!\times\binom{n+m-1}{m-1}=\frac{(n+m-1)!}{(m-1)!}\ . $$ This proves the more general conjecture in Dan Piponi's comment.