Let's step back and ask what the determinant is and why it is useful. Then we'll get to this particular property.
Permit me to introduce a new product of vectors, called the wedge product. If $a, b$ are vectors, then their wedge product is $a \wedge b$. The result is not a vector, but instead, we interpret it geometrically as an oriented plane--exactly the plane perpendicular to $a \times b$, as a matter of fact, but we can continue wedging--for instance, forming $a \wedge b \wedge c$, which represents a volume.
Of course, in 3d space, there is only one unit volume (you could call it $\hat x \wedge \hat y \wedge \hat z$, but let me call it $i$ for short. All other volumes are scalar multiples of this volume, at least in terms of magnitude and orientation (you may be thinking "orientation?" I submit that a coordinate system following the right-hand rule is oppositedly oriented to one following a left-hand rule, and the unit volumes formed by wedging their unit vectors are oppositely oriented). It is for this reason that the volume elements are often called pseudoscalars, because they are only different from scalars by how they have orientations where scalars do not.
How does this relate to matrices? Well, matrices are used to represent linear operators on vectors, but these operators can act on wedge products of vectors or pseudoscalars too. The "matrices" you use to represent these extensions of the original operator are different from the original matrix. We define the relationships by a simple rule. If $\underline T$ is a linear operator on a vector, then we define $\underline T(a \wedge b) = \underline T(a) \wedge \underline T(b)$, and so on (note that the cross product does not have this nice property except under rotations).
But again, we said that there is only one unit pseudoscalar ($i$), so the action of a linear operator on $i$ must be some scalar multiple of $i$. That is, if $\alpha$ is a scalar,
$$\underline T(i) = \alpha i$$
We define this number $\alpha$ to be the determinant, telling us geometrically how the unit volume is shrunk or dilated (or changes orientation) under the action of a linear operator.
You can find the determinant of a linear operator by writing down the matrix representation and wedging the vectors that appear there. This is perfectly well-founded. You just need to know that $a \wedge b = - b \wedge a$--vectors anticommute under the wedge--and that the wedge is associative. Knowing these properties makes it possible to do computations with the wedge.
(If you're puzzled why wedging vectors gives the determinant, feel free to ask and I'll clarify this point.)
So let's take three vectors $f, g, h$ that appear in the matrix representation of a linear operator and wedge them to find the determinant. Let $f = f^x e_x + f^y e_y + f^z e_z$ and so on. We can then write out the following:
$$\begin{align*}
f \wedge g \wedge h = f^x e_x \wedge (g \wedge h) + f^y e_y \wedge (g \wedge h) + f^z e_z \wedge (g \wedge h)\end{align*}$$
I've expanded the wedge product through linearity (the distributive property, which I failed to mention earlier but is also valid for the wedge). This is the foundation for the technique you'e come across, which is called Laplace expansion, or cofactor expansion, or expansion by minors. The antisymmetry of the wedge means that we can write $f^x e_x \wedge (g \wedge h)$ as
$$f^x e_x \wedge (g \wedge h) = f^x e_x \wedge (g^y e_y + g^z e_z) \wedge (h^y e_y + h^z e_z)$$
Why? Because $e_x \wedge e_x = 0$ always, so if any $e_x$ appeared in $g$ or $h$, they would be irrelevant. We can ignore them, and instead, we find the "determinant" of a linear operator on the $yz$-plane. Thus, the method of expansion by minors follows by recursively expanding a single vector through linearity to reduce the subsequent wedge products to a wedge product (a "determinant") you might already know. The method has its roots in the relationship between the determinant and volumes, how all the terms that make up the determinant must contain 1 and only 1 component of some vector from each of the coordinate directions--no more, no less.
Best Answer
As you pointed out, the integrand is invariant under the action of the orthogonal group. Thus it suffices to integrate over diagonal matrices, and multiply the result by the volume of the orthogonal group.
Suppose $X=diag(x_1,\ldots,x_k)$. Then the integral becomes $$ \int_{\mathbb R} \prod_{i=1}^k\left(1+\frac{g}{a}x_i^2\right)^{-d/2}e^{-ax_i^2/2}\ dx_1\cdots dx_k=\left[\int_{\mathbb R}\left(1+\frac{g}{a}x^2\right)^{-d/2}e^{-ax^2/2}\ dx\right]^k. $$ The innermost integral may be expressed in terms of moments of the Gaussian distribution.