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 therem 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 \\ | \\ X \\ | \end{array} \begin{array}{c} | \\ X \\ | \\ X \\ | \end{array} \right]= \begin{array}{} | & & | \\ - & -- & - \\ & & \\ - & -- & - \\ | & & | \end{array} $$ or with indices $$ \mathbb{E}[X_{ij}X_{k\ell}]=\delta_{ik}\delta_{j\ell} $$ 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)$. This gives 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.