Generating function of $SO(N)$ random matrix

For the $n=3$ case, it is best to use quaternionic methods to express $SO(3)$ as a quotient of $S^3$. Explicitly, for $a\in S^3$ one can check that the matrix $$ A = \left[\begin{matrix} a_1^2-a_2^2-a_3^2+a_4^2 & 2a_1a_2-2a_3a_4 & 2a_1a_3+2a_2a_4 \\ 2a_1a_2+2a_3a_4 & -a_1^2+a_2^2-a_3^2+a_4^2 & -2a_1a_4+2a_2a_3 \\ 2a_1a_3-2a_2a_4 & 2a_1a_4+2a_2a_3 & -a_1^2-a_2^2+a_3^2+a_4^2 \end{matrix}\right] $$ lies in $SO(3)$. (This is based on an identification $\mathbb{H}\simeq\mathbb{R^4}$ with $1\in\mathbb{H}$ corresponding to $(0,0,0,1)\in\mathbb{R}^4$.) This construction comes from a group homomorphism, so the Haar measure on $SO(3)$ corresponds to the natural round measure on $S^3$, suitably rescaled. Thus, we just want to integrate the function $\exp(f)$ over $S^3$, where $$ f(a) = (a_1^2-a_2^2-a_3^2+a_4^2)\lambda_1 + (-a_1^2+a_2^2-a_3^2+a_4^2)\lambda_2 + (-a_1^2-a_2^2+a_3^2+a_4^2)\lambda_3. $$ Using $\sum_ia_i^2=1$, we get $f(a)=\sum_i\lambda_i-2g(a)$, where $$ g(a) = (a_2^2+a_3^2)\lambda_1 + (a_1^2+a_3^2)\lambda_2 + (a_1^2+a_2^2)\lambda_3.$$ This is at least visibly invariant under permutations of $\lambda$. Next, because $a_4$ does not appear in $g(a)$, it is natural to use a stereographic parametrisation $\sigma\colon\mathbb{R}^3\to S^3$, as follows: $$ \sigma(x_1,x_2,x_3) = \frac{(2x_1,2x_2,2x_3,\|x\|^2-1)}{\|x\|^2+1} $$ The Jacobian determinant of $\sigma$ turns out to be $8/(\|x\|^2+1)^3$. This leads to an expression for $Z_3$ as an integral over $\mathbb{R}^3$ that retains all the natural symmetries, but it does not seem easy to evaluate either numerically or symbolically.

Singular value decomposition: The real $N\times N$ matrix $J$ has singular values $\sigma_1,\sigma_2\ldots\sigma_N\geq 0$ in the decomposition $J=U\,{\rm diag}\,(\sigma_1,\sigma_2\ldots\sigma_N)V$ with $U,V\in{\rm O}(N)$. We need to restrict $U,V$ to ${\rm SO}(N)$, which we can do by allowing for one of the singular values to become negative: $$J=U\,{\rm diag}\,(\lambda_1,\lambda_2\ldots\lambda_N)V,\;\;U,V\in{\rm SO}(N),$$ $$\lambda_n=\sigma_n\geq 0,\;\;{\rm for}\;\;n=1,2,\ldots N-1,$$ $$\lambda_N=\begin{cases} \sigma_N\;{\rm if}\;{\rm det}\,M\geq 0,\\ -\sigma_N\;{\rm if}\;{\rm det}\,M<0. \end{cases}$$ Note that this need to consider the $\lambda_n$'s implies that the integral $Z_N$ is not only a function of the invariants ${\rm Tr}((J^TJ)^n)$.

Without loss of generality we may now take a diagonal $J={\rm diag}\,(\lambda_1,\lambda_2,\ldots \lambda_N)$ and then we need the integral $$Z_N=\int_{{\rm O}(N)} dM \exp\left(\sum_{n=1}^N \lambda_n M_{nn}\right).$$
Large-$N$ limit: I do not think there is a closed form expression for arbitrary $N$, but for $N\gg 1$ the diagonal matrix elements of $M$ are independent Gaussians (mean zero, variance $1/N$), resulting in $$Z_N\rightarrow \exp\left(\frac{1}{2N}\sum_{n=1}^N \lambda_n^2\right),\;\;N\gg 1.$$ This result for $Z_N$ in the large-$N$ limit depends only on the singular values (because $\lambda_n$ enters only quadratically), but for small-$N$ this is no longer true.

Small-$N$ results: For small $N$ we can use the parametrization of the orthogonal group given in Appendix B of arXiv:1405.3115.

• $N=2$ $$Z_2=\frac{1}{\pi}\int_0^{\pi} \exp\bigl[(\lambda_1+\lambda_2)\cos\theta\bigr]\,d\theta=I_0(\lambda_1+\lambda_2)$$ Note that in terms of the singular values considered in the OP, this would read $Z_2=I_0(\sigma_1\pm\sigma_2)$, where the $\pm$ sign is the sign of ${\rm det}\,J$.

• $N=3$ $$Z_3=\frac{1}{8\pi^2}\int_0^{2\pi} d\alpha\int_0^{2\pi} d\alpha'\int_0^\pi \sin\theta\, d\theta$$ $$\qquad\qquad\times\exp\left[\cos \alpha \cos \alpha' (\lambda_1+\lambda_2 \cos \theta)-\sin \alpha \sin \alpha' (\lambda_1 \cos \theta+\lambda_2)+\lambda_3\cos \theta\right]$$ The two integrals over $\alpha,\alpha'$ can be carried out (by rewriting them as integrals over $\alpha\pm\alpha'$), with the result $$Z_3=\frac{1}{2}\int_{-1}^1 I_0\left[\tfrac{1}{2}(\lambda_1-\lambda_2) (1-x)\right] I_0\left[\tfrac{1}{2} (\lambda_1+\lambda_2)(1+x)\right]\,e^{\lambda_3 x}\,dx.$$ As shown by Paul Enta, this integral can be written as the inverse Laplace transform ${\cal L}^{-1}[F(s)](t)$ for $t=1$ of the function $$F_3(s)=\frac{1}{\sqrt{(s-\lambda_1-\lambda_2-\lambda_3) (s+\lambda_1+\lambda_2-\lambda_3) (s+\lambda_1-\lambda_2+\lambda_3) (s-\lambda_1+\lambda_2+\lambda_3)}},$$ to demonstrate that it is invariant under permutation of the $\lambda_n$'s.
_{Notice also that $F_3(s)$ is invariant under a sign change of two of the $\lambda_n$'s, but not under a sign change of one single $\lambda_n$. This is consistent with the definition of the $\lambda_n$'s in terms of the singular values, given above.}

A closed-form expression of either the Bessel-function integral or the inverse Laplace transform seems not forthcoming, except for some special cases (see the comments by Adam).

Average over the full orthogonal group: If we integrate over the full group ${\rm O}(N)$, instead of only over ${\rm SO}(N)$, we should combine the results for $\pm\lambda_N$, $$Z_{{\rm O}(N)}(\lambda_1,\lambda_2,\ldots \lambda_{N-1},\lambda_N)= \tfrac{1}{2}Z_{{\rm SO}\,(N)}(\lambda_1,\lambda_2,\ldots \lambda_{N-1},\lambda_N)+ \tfrac{1}{2}Z_{{\rm SO}(N)}(\lambda_1,\lambda_2,\ldots \lambda_{N-1},-\lambda_N).$$ There is now no need to distinguish the $\lambda_n$'s from the singular values $\sigma_n$. In particular, for $N=2$ one has $$Z_{{\rm O}(2)}=\tfrac{1}{2}I_0(\sigma_1+\sigma_2)+\tfrac{1}{2}I_0(\sigma_1-\sigma_2).$$