What happens is similar to when, in the 1D case, the variance is zero: the variable becomes deterministic, so you would need to use "unorthodox tricks" (Dirac delta) to characterize the random variable with a Probability Density Function (PDF)*.
What becomes deterministic in the $N$-D case is a set of linear combinations of the random variables, and the number of these deterministic combinations (as in: dimension of the corresponding linear space) is determined by the rank of the matrix (more precisely, its deficiency).
*But refer to @Michael Hardy's answer for the proper characterizations which require neither the existence of a PDF nor the use of a Dirac delta.
I believe this is simpler to intuitively explain if we restrict to the case of 2D/bivariate jointly Gaussian distributions. So, consider two random variables $X$ and $Y$.
To begin with, it is important to remark that we are only interested in the case where $X$ and $Y$ are jointly Gaussian, which means that any linear combination of the variables is Gaussian (1D) or a deterministic constant; the latter can be considered as a special case of a ("generalized") Gaussian where the variance tends to zero (this is actually relevant for the problem at hand).
Note that $(X,Y)$ being jointly Gaussian is a stronger condition than just having both $X$ Gaussian and $Y$ Gaussian separately (i.e. the marginals are Gaussians); we have to assume this restriction because otherwise several unintuitive issues would arise (e.g. the "uncorrelated and Gaussian implies independence" rule is broken).
Now, under this jointly Gaussian assumption, the set of possible PDFs is constrained; in particular, for the 2D case, and assuming on top of it that the variables are not linearly dependant, the PDF has to have this structure:
$$
f_{XY}(x,y) =
\frac{
\text{exp}
\left(
-\frac{1}{2(1-\rho^2)}\left[
\left(\frac{x-\mu_X}{\sigma_X}\right)^2 -
2\rho\left(\frac{x-\mu_X}{\sigma_X}\right)\left(\frac{y-\mu_Y}{\sigma_Y}\right) +
\left(\frac{y-\mu_Y}{\sigma_Y}\right)^2
\right]
\right)
}{2 \pi \sigma_X \sigma_Y \sqrt{1-\rho^2}}
$$
You may identify the structure from the general case ($N$-dimensional jointly Gaussian):
$$
f_{\mathbf X}(x_1,\ldots,x_N) = \frac{\exp\left(-\frac 1 2 ({\mathbf x}-{\boldsymbol\mu})^\mathrm{T}{\boldsymbol\Sigma}^{-1}({\mathbf x}-{\boldsymbol\mu})\right)}{\sqrt{(2\pi)^N |\boldsymbol\Sigma|}}
$$
Where the matrices that parametrize the 2D distribution can be written as:
$$
\boldsymbol\mu = \begin{pmatrix} \mu_X \\ \mu_Y \end{pmatrix}, \quad
\boldsymbol\Sigma = \begin{pmatrix} \sigma_X^2 & \rho \sigma_X \sigma_Y \\
\rho \sigma_X \sigma_Y & \sigma_Y^2 \end{pmatrix}
$$
And with this, it can be proven that the two PDF expressions are indeed equivalent.
Now let's focus on the parameter $\rho\,$: this is the so-called Pearson correlation coefficient and indicates a degree of correlation between two random variables, it is defined as
$$\rho =
\rho_{X,Y}= \frac{\operatorname{cov}(X,Y)}{\sigma_X \sigma_Y}
= \frac{\mathbb{E}[(X-\mu_X)(Y-\mu_Y)]}{\sigma_X\sigma_Y}
$$
$\rho$ ranges from $-1$ to $+1$.
One can see from the first PDF expression that the extreme values $\pm 1$ are problematic: this actually corresponds to the case when the covariance matrix $\boldsymbol\Sigma$ is singular. This can be verified with the determinant:
$$|\boldsymbol\Sigma| = \sigma_X^2\sigma_Y^2-\rho^2\sigma_X^2\sigma_Y^2=\sigma_X^2\sigma_Y^2(1-\rho^2)=0$$
To see what is actually happening, let's start with the opposite case: $\rho=0$, this corresponds to uncorrelated variables; they are also independent, as they are jointly Gaussian:
$$\rho=0\implies f_{XY}(x,y)
=\frac{
\text{exp}
\left(
-\frac{1}{2}\left[
\left(\frac{x-\mu_X}{\sigma_X}\right)^2 +
\left(\frac{y-\mu_Y}{\sigma_Y}\right)^2
\right]
\right)
}{2 \pi \sigma_X \sigma_Y}
= f_{X}(x)\cdot f_Y(y)
$$
For $\rho=0$, we have independence, which can be observed if we plot the PDF (3D representation): we get a circular pattern if $\sigma_X=\sigma_Y$, elliptic otherwise. For the sake of example let's assume that $\sigma_X=\sigma_Y$ (this also implies slope $\pm 1$ as $\rho\to \pm 1$) and $\mu_X=\mu_Y=0$ (symmetry center at $(0,0)$).
Here we see how the changes in $\rho$ affect the resulting distribution (from Wikipedia):
If $\rho\neq 0$, there is a degree of linear dependence between $X$ and $Y$ that is made perfect $X=\pm Y$ at the extreme values of $\rho$ (deterministic dependence). Note that the sign of $\rho$ establishes the sign of the slope.
Now we can see what the problem is: $f_{XY}(x,y)$ is 2D density that is supposed to be integrated by a double integral in order to produce a probability. And also note that the total probability over $\mathbb R^2$ must be 1:
$$
\iint\limits_{\mathbb R^2} f_{XY}(x,y)\,dx\,dy = 1
$$
But in the problematic extreme cases, the only region that is nonzero is a straight line, and thus it has zero area (zero measure). So, to add up to 1, the density at those line points must be "infinite"! This is the reason why we cannot define $f_{XY}(x,y)$ as an actual function. On the other hand, note that it would still make sense to define a 1D density function on that line, according to a Gaussian PDF.
To overcome this problem, we can make use of a "generalized" function (distribution), the Dirac delta: this provides a method to define a density that, in principle, should be infinite, but is defined in such a way that it makes sense under integration. The effective result is that the region of integration is reduced from a 2D region ($\mathbb{R}^2$ plane) to a 1D region (line, $y=\pm x$ for $\rho=\pm 1$, respectively).
Therefore, we can now define the PDF, e.g. for $\rho=1$:
$$f_{XY}(x,y) = \delta(y-x) \cdot f_X(x)$$
where $f_X(x)$ is the 1D Gaussian PDF.
With that, we defined a 2D PDF from a density that is in reality 1D:
$$
f_X(x) =
\int_{-\infty}^{+\infty}
f_{XY}(x,y)\,dy =
f_X(x)\cdot
\int_{-\infty}^{+\infty}
\delta(y-x)\,dy
= f_X(x)\cdot 1
$$
And therefore:
$$
\iint\limits_{\mathbb R^2}f_{XY}(x,y) \,dxdy =
\iint\limits_{\mathbb R^2}\delta(y-x) \cdot f_X(x) \,dx\,dy =
\int\limits_{\text{line } y=x}f_X(x) \,dx = 1
$$
Now, following the argument in the first paragraph: for $\rho=1$ there is a set of non-trivial linear combinations that are deterministic: $Z = a\cdot X - a\cdot Y$ for all $a\in \mathbb{R}$.
In this case, $Z=0$ but it could be any deterministic set of constants (that depends on the parameters $\boldsymbol\mu,\, \boldsymbol\Sigma$ and on $a$).
In the general $N$-D case, the presence of such combinations is a symptom of the described problem: the density is actually "infinite" and the PDF doesn't exist as a proper function; it requires a Dirac delta to be defined. Again, this will convert the probability integral "as if" the number of variables/dimensions was smaller. And, in particular, the number of dimensions that has to be reduced is precisely the deficiency of the covariance matrix ($N-\text{rank}(\Sigma)$, i.e. number of dimensions away from being full rank).
