How to generate random samples of Gaussian distribution directly in the frequency domain?
You can, but... you'll need to keep symmetry if your original time-domain signal is real-valued.

If a signal $x$ is real-valued, then its DFT $X$ will exhibit complex-conjugate symmetry:
$$
X[k] = X^*[N-k].
$$

So you can generate $N$ Gaussian pseudo-random noise samples, $g[n]$, and place them in the frequency domain noise vector, $\epsilon$ as:
$$
\epsilon[k] =  g[k] + j g[k+N/2] 
$$
for  $k \in \{ 0, 1, \ldots, N/2-1\}$ and
$$
\epsilon[k] = \epsilon^*[N-k] 
$$
for $k \in \{ N/2, N/2+1, \ldots, N-1 \}$ where $\epsilon^*$ is the complex conjugate of $\epsilon$ and is equal to $\Re[{\epsilon}] - \Im[{\epsilon}]$ (i.e. the same real part and the negative of the imaginary part).