Q:

A general formula for Gaussian integrals over matrix elements

The question I have is quite specific.  So in the hope that this post might help others in the future, my problem boils down to solving the following integral:
$$I_\tau=\int \prod_{i, j=1}^{N} d J_{i j} \exp  \left\{-\frac{N}{2(1-\tau^2)}\left( \sum_{i, j, k} J_{k i} A_{i j} J_{k j}-\tau\sum_{ij}J_{ij}J_{ji}\right)\right\}$$
Assuming we can compute the determinant of $A$, is there a general identity to compute $I_\tau$? 

Simplifying their argument, in the following paper [1], the authors compute the following Gaussian integral over all $N\times N$ real asymmetric matrices $J_{ij}$: 
$$I=\int \prod_{i, j=1}^{N}d J_{i j} \exp \left\{-\frac{N}{2} \sum_{i, j, k} J_{k i} A_{i j} J_{k j}\right\}$$
Where $A_{i j}=\frac{z_{i}^{*} z_{j}}{N}+\frac{z_{j}^{*} z_{i}}{N}+\delta_{i j}$ and $z_i\in \mathbb{C}$.
If I define the vector $\mathbf{x}$ such that:
$$\mathbf{x}:=\left(\begin{array}{c}
J_{11} \\
J_{12} \\
\vdots \\
J_{1 n} \\
J_{21} \\
J_{22} \\
\vdots \\
J_{n n}
\end{array}\right) \in \mathbb{R}^{n^{2}}$$
$$\implies x_{N(i-1))+j}=J_{ij}$$
Then we can represent $I$ as the following integral:
$$I=\int \left(\prod_{i, j=1}^{N} \mathrm{d}x_{N(i-1))+j} \right)\exp \left(-\frac{N}{2} \mathbf{x}^T\Sigma \mathbf{x}\right)$$
Where $\Sigma$ simply is:
$$\Sigma=\left(\begin{array}{cccc}
\mathbf{A} & 0 & \cdots & 0 \\
0 & \mathbf{A} & \cdots & 0 \\
\vdots & \vdots & \ddots & \vdots \\
0 & 0 & \cdots & \mathbf{A}
\end{array}\right)$$
Thus $I\propto \frac{1}{\sqrt{|\Sigma|}}=|A|^\frac{-N}{2}$.
Defining $r\equiv\frac{1}{N}\sum_i |z_i|^2$, the authors observe that to order $1/\sqrt{N}$, the matrix $A$ has all but two of its eigenvalues equal to one.  The two exceptions are the eigenvectors $z_i$ and $z^*_i$ , both with eigenvalue $(1+r)$. This observation allows them to compute the determinant of $\Sigma$ and thus the Gaussian integral $I$.

Now, considering the integral $I_\tau$ using the same method:
$$I_\tau=\int \left(\prod_{i j} \mathrm{d} J_{i j}\right)\exp  \left\{-\frac{N}{2(1-\tau^2)}\left( \sum_{i, j, k} J_{k i} A_{i j} J_{k j}-\tau\sum_{ij}J_{ij}J_{ji}\right)\right\}$$
$$I_\tau=\int\left(\prod_{i j} \mathrm{d} J_{i j}\right) \exp \left\{-\frac{N}{2(1-\tau^2)} \sum_{i j, k, l} J_{i j}\left(A_{j l} \delta_{i k}-\tau\delta_{i l} \delta_{j k}\right) J_{k l}\right\}$$
$$\implies I_\tau=\int\left(\prod_{i, j=1}^{N} \mathrm{d} x_{N(i-1))+j}\right) \exp \left(-\frac{N}{2} \mathbf{x}^{T} \Sigma \mathbf{x}\right)$$
Where this time:
$$\Sigma_{n(i-1)+j,n(k-1)+l}=A_{j l} \delta_{i k}-\tau\delta_{i l} \delta_{j k}  \quad \forall i, j, k, l \in \mathbb{N} \cap[1, N]$$
For large $N$, how can we compute the determinant of $\Sigma$ this time?  It seems to be a non-trivial question. 
Any remark or comment is always appreciated, thank you. 

[1]: Rajan, K., & Abbott, L. F. (2006).  Eigenvalue spectra of random matrices for neural networks.  Physical review letters, 97(18), 188104. 
Link: https://www.researchgate.net/profile/Kanaka_Rajan/publication/6643146_Eigenvalue_Spectra_of_Random_Matrices_for_Neural_Networks/links/0912f50f7347f4863b000000.pdf

A:

Only the symmetric part of $A$ contributes to the integrand, so we may assume $A$ is symmetric and diagonalize it as $A=O\Lambda O^T$ with $O$ orthogonal and $\Lambda={\rm diag}\,(\lambda_1,\lambda_2,\ldots\lambda_N)$. The orthogonal matrix may be removed by a change of variables, so without loss of generality we can take a diagonal $A$.  So we seek the integral
$$I=\int \prod_{i, j=1}^{N} d J_{i j} \exp  \left\{-\frac{N}{2(1-\tau^2)}\left( \sum_{i, j} \lambda_j J_{ij}^2-\tau\sum_{ij}J_{ij}J_{ji}\right)\right\}.$$
I now proceed similarly to here.  Decompose the sum over $i,j$ as
$$\sum_{ij}\lambda_j J_{ij}^2-\tau\sum_{ij}J_{ij}J_{ji}=$$
$$\qquad\qquad=\sum_{i}(\lambda_i-\tau)J_{ii}^2+\sum_{i<j}\left(\lambda_j J_{ij}^2+\lambda_i J_{ji}^2-2\tau J_{ij}J_{ji}\right)$$
$$\qquad\qquad\equiv\sum_{i} A_i+\sum_{i<j} B_{ij}.$$
Then perform the Gaussian integrals separately for each term in the sum, 
$$I=\left(\prod_{i=1}^N\int e^{-\beta A_i}dJ_{ii}\right)\left(\prod_{i<j=1}^{N}\int\int e^{-\beta B_{ij}}dJ_{ij}dJ_{ji}\right)$$
$$\qquad\qquad=(\pi/\beta)^{N^2/2}\left(\prod_{i=1}^N(\lambda_i-\tau)^{-1/2}\right)\left(\prod_{i<j=1}^N(\lambda_i\lambda_j-\tau^2)^{-1/2}\right),$$
where I have defined $\beta=\frac{1}{2}N(1-\tau^2)^{-1}$, and assumed that $\beta>0$, $\lambda_i>\tau$ for all $i$. 
So knowledge of only the determinant of $A$ is not sufficient to evaluate the integral, you need to know the individual eigenvalues -- not just their product.  In particular, if $\tau$ happens to be close to one particular eigenvalue of $A$, then that eigenvalue will dominate the integral.