Your question is an interesting project for you to research it on your own. Well, maybe with a little help from your friend and the SE community. And, as your question goes, start with generating samples from (1) a uniform distribution, and (2) a Gaussian distribution. I will accompany you in your first steps.
First, I generate waveforms of 1,048,576 samples for each type of noise. Here in the pictures the ranges of 512 samples long are shown.

Uniform noise

Gaussian noise
The pictures look quite similar.
Then I compute FFTs of waveforms for each of type of noise and plot the distributions of values of the noise waveform (blue graph) and of the corresponding FFT (orange graph).

Uniform noise. Distribution of sample values in blue, distribution of FFT values in orange

Gaussian noise. Distribution of sample values in blue, distribution of FFT values in orange
I do not quite understand what you mean talking about "a flat frequency spectrum". Both signals are realizations of stochastic processes. But here they are, and the trends in the graphs of sample value distributions clearly betray their origin. At the same time, distributions of FFT values (real parts of FFT values) for both signals have a distinct look of Gaussian distribution. To prepare yourself to follow the explanation of this behavior, I recommend you to study one of the proofs of the central limit theorem of probability theory, the proof, which uses Fourier transform. With this knowledge you may even find the explanation on your own. If not, still it is a useful thing to know when following the explanations in the literature.
Notice also, in the graphs of the distribution of FFT values, the central peaks, which are two times higher then the interpolated distribution at zero.
UPDATE
While modeling white noise in computations, under certain conditions you are free to select any method that will generate random values for samples in your simulation. The generated values can follow uniform, Gaussian, any other distributions with a zero mean value. Even a coin toss random generator, randomly producing +ampl/-ampl with a 50/50 chance, will fit. Of these "certain conditions", one is just mentioned -- a zero mean value. Another requirement is a high sampling frequency, the higher the better. Within any sample interval between adjacent signal samples you must generate a great many noise samples, densely filling the time interval. This, admittedly nonobvious, procedure is justified by the essence of "white noise signal": it is a stochastic process, not a classical function. In an ordinary language, it means that we cannot know a noise value at any given moment, but noise values averaged over an arbitrarily small time intervals can be "measured" in experiment, or "simulated" in computations following a known formula, this "known formula" following from the central limit theorem of probability theory.
With this approach, the signal sampling (and, consequently, Nyquist) frequency is much lower than the noise sampling frequency. To avoid aliases in the noise frequency domain, you must bandlimit the noise. In this process, the noise is averaged over the time interval between the adjacent signal samples in a manner dependent of the filter used. By force of the central limit theorem, the averaged noise values tend to approach a Gaussian distribution -- the more samples used in averaging process, the closer the distribution to Gaussian. But this fact suggest the idea to directly generate samples from a Gaussian distribution and avoid oversampling at frequency much higher than the signal sampling frequency.
The method generating a random variable with a Gaussian distribution is readily available in matlab, numpy and scipy (I believe). In the general-purpose computer languages you can easily implement the Marsaglia polar method or the Box–Muller transform.
Summing up, the "It turns out that..." statement from you question is the corollary of the central limit theorem. You can see it either following the central theorem proof or even better generating a number of i.i.d. random variables and computing the distribution of their sum, repeating the process with an increased number of random variables and examining the results.
- The generation of noise via (1) never suffers from aliasing in principle, but the computed spectra can be aliased if an inadequate Nyquist frequency is used for computation. With the Nyquist frequency fixed by the signal sampling frequency, you must adjust the noise sampling frequency.
- Aliasing manifests itself via a phenomenon called Spectral Leakage.
- The standard deviations resulting in the signal values exceeding the maximum values of the variable types used in the program do chop the tail of the distribution. Using float's rather than UINT16's is a safe option to fight this chopping.
NUMERICAL EXPERIMENT: NOISE GENERATION
Start with a silent (zero-valued) signal sampled into a data array of 65536 samples. First, add a uniform-distributed noise, one noise sample per one signal sample. The first 256 samples of realization (pure noise) are shown here:

Compare the distribution of sample values against a Gaussian distribution of equal power:

The FFT of the realization:

The PSD of the realization:

The zoomed in PSD of the realization 256 samples long:

Back to the original zero-valued pure signal, and add a uniform-distributed noise again, but this time insert two noise samples per signal sample. To produce a realization of the noisy signal sampled at the original pure signal Nyquist frequency, we sum the adjacent pairs of samples. To keep the noise power constant, we scale down the noise amplitude by the square root of oversampling (=2). The first 256 samples of the new realization are shown here:

Compare the distribution of sample values against a Gaussian distribution of equal power:

With just two-fold oversampling, the distribution is already close to the Gaussian distribution, as it is seen in the difference plot of the two distributions:

The zoomed in PSD of the realization 256 samples long:

Not much different for the plots of purely uniform noise.
Finally, calculate with the oversampling of 4096 noise samples per signal sample:
The realization after averaging (first 256 samples):

The difference of the two distributions, uniform-folded and Gaussian:

The PSD of the realization:

UPDATE on AWGN:
The phrase "but the computed spectra can be aliased if an inadequate Nyquist frequency is used for computation" may lead to confusion. I'll try to explain it without explicitly referring aliasing (the aliasing thing still implicitly leaks through the references, but these are presumably more trustworthy than my explanations).
See it the other way round: when, in your attempt to generate noise, you generate one random uniform-distributed variable per sound sample (the sampling frequency for generated noise is 44.1KHz), the PSD is also random and very uneven. To approximate the white noise's PSD behavior, you can resort to N-oversampling, increasing the sampling and Nyquist frequency to N·44.1KHz: you spread each value of sound sampled at 44.1KHz among adjacent N samples of a new (N·sample_count)-sized vector (array) of noised samples and add uniform noise.
std::random_device rd;
std::mt19937 gen(rd());
std::uniform_real_distribution<double> unirand(-1.0, 1.0);
for(int isample=0; isample < vecSig.size()/*44.1-samples_count*/; ++isample)
for(int ix=0; ix < N; ++ix)
vecOversampled[isample*N + ix] = vecSig.at(isample) + unirand(gen);
The PSD of this signal averaged over N adjacent frequency bins is a way more even. Seeing that the frequencies higher than 44.1KHz are sort of computation crutches, we filter out these frequencies from the signal and arrive at the equivalent of the Gaussian-distributed random variables added to the original 44.1KHz samples, because the central limit theorem holds for the averaged sample values. Alas, the PSD of the filtered signal is ragged anew, for the PSD computation is a non-linear operation, and the PSD plot can be smoothed only with enlarged frequency bin sizes. But this time we know that this raggedness is due to windowing function and we have a leaking in process rather than leaking out.
This PSD raggedness is the result of sampling of the white noise "signal".
Summing up: to spare memory and CPU time, you can do without oversampling, generate a Gaussian-distributed noise from the beginning
std::random_device rd;
std::mt19937 gen(rd());
std::normal_distribution<double> normalrand(-1.0/std::sqrt(M_PI), 1.0/std::sqrt(M_PI));
for(int isample=0; isample < vecSig.size()/*44.1-samples_count*/; ++isample)
vecSigPlusNoise[isample] = vecSig.at(isample) + normalrand(gen);
and arrive at the same results as with the oversampled uniform-distributed noise. You cannot do without oversampling when modeling the noise with uniform-distributed random variables -- without delving into the depths of probability theory, I refer you to the Additive White Noise definition (https://en.wikipedia.org/wiki/White_noise)
In discrete time, white noise is a discrete signal whose samples are regarded as a sequence of serially uncorrelated random variables with zero mean and finite variance; a single realization of white noise is a random shock. Depending on the context, one may also require that the samples be independent and have identical probability distribution (in other words independent and identically distributed random variables are the simplest representation of white noise). In particular, if each sample has a normal distribution with zero mean, the signal is said to be additive white Gaussian noise. [my emphasis]
whereas each sample in a non-oversampled uniform-distributed noise model has a uniform distribution.
I presume that at this stage you need an AWGN simulation for your project.
It may be just an intuition in the requirement that the "ideal" noise model must exhibit similar patterns at increasingly small scales, a so-called self-similarity feature. You see that if you start with a uniform-distributed or any other non-Gaussian-distributed noise at some fine frequency resolution and "zoom out" this noise to a larger, coarser scale, you arrive at the noise tending to approach a Gaussian-distributed noise -- whereas the Gaussian-distributed noise reproduces the Gaussian distribution of sample values at any zoomed resolution (with a properly adjusted standard deviation parameter). This self-similarity feature is often used for the analysis and computation in noise-related development projects. As soon the white noise model is applicable, you can safely assume it is AWGN in discrete time and use the self-similarity feature for analysis.
In addition to the self-similarity of noise value distributions among zoomed frequency resolutions, another argument for Gaussian distribution versus any other distribution can be derived from a process of direct construction of Fourier spectrum to satisfy a "flat PSD" requirement.
The straightforward function of frequency which gives the most flat, indeed constant, PSD, is a constant, i.e., independence of frequency. Unfortunately, the signal, of which the Fourier spectrum is independent of frequency, is a delta function, and the delta function is everything but the noise.
To endow the constant Fourier spectrum with noise-like traits, one can follow your recipe for producing the noise in time domain: each sample of frequency spectrum is made a random uniform-distributed variable. But, as we have seen, the Fourier spectrum constructed in this way gives the time-domain realization with Gaussian-distributed sample values!
We have discussed two possibilities to simulate white noise. In one approach, we start with the uniform-distributed samples in a time domain and have the Gaussian-distributed samples in a frequency domain. In another approach, we start with the uniform-distributed samples in a frequency domain and have the Gaussian-distributed samples in a time domain. Interestingly, this means that there exist realizations with Gaussian-distributed samples in the frequency domain for which the procedure of inverse Fourier transform gives uniform-distributed samples in the time domain. But the share of these realization in the signal space is negligible (a set of measure zero), and for an overwhelming majority of signal realizations, starting with the Gaussian-distributed sample values in the time (frequency) domain, one obtains the Gaussian distributed values in the frequency (time) domain.
The Gaussian distribution appears to be the most natural distribution for generating white noise realizations. In the signal space, the subspace of non-Gaussian-distributed white noise realizations has measure zero.