The purpose of this answer is simply to expand on the answer by @Tim.
Suppose the likelihood of the parameters given the sample can be expressed as
\begin{equation}
p(X|\theta) = \prod_{i=1}^n \textsf{N}(x_i|\mu,\sigma^2) ,
\end{equation}
where $X = (x_1, \ldots, x_n)$ is the sample and $\theta = (\mu,\sigma^2)$ are the parameters. Then in general the likelihood of model $j$ (i.e., class $j$) can be expressed as
\begin{equation}
p(X|C_j) = \int p(X|\theta)\,p(\theta|C_j)\,d\theta ,
\end{equation}
where $p(\theta|C_j)$ is the distribution of $\theta$ given model $j$.
This general approach can be specialized to the current case as follows. Let
\begin{equation}
p(\theta|C_j) = \delta(\mu-\mu_j)\,\delta(\sigma^2 - \sigma_j^2) ,
\end{equation}
where $\delta(x)$ is the Dirac delta "function." In effect, this distribution puts a point mass at $\theta_j = (\mu_j,\sigma_j^2)$. The two salient properties of the Dirac delta function are $\int \delta(x-x_0)\,dx = 1$ and $\int f(x)\,\delta(x-x_0)\,dx = f(x_0)$.
With this point-mass distribution, we can compute the desired expression:
\begin{equation}
\begin{split}
p(X|C_j) &= \iint p(X|\mu,\sigma^2)\,\delta(\mu-\mu_j)\,\delta(\sigma^2-\sigma_j^2)\,d\mu\,d\sigma^2 \\
&= p(X|\theta_j) \\
&= \prod_{i=1}^n \textsf{N}(x_i|\mu_j,\sigma_j^2) .
\end{split}
\end{equation}