Why are the observation features of an HMM-based recognition/synthesis system modeled by a Gaussian distribution?
For observations people usually use a mixture of gaussians, not a simple gaussian. They have few advantages - EM algorithm is fast to converge and GMM approximates wide variety of distributions pretty well. Probability of the model can be computed efficiently with gaussian selection. Last, GMMs are easy to cluster for context-dependent tree for phonetic context, which is important for LVCSR.

There are disadvantages as well, for example, GMM estimation is not very robust to the noise in features, also it does not account for some hidden parameters, for example "speaker" identity is not properly modeled with GMM. 

There are extensions developed here, for example, subspace GMM is proposed to model speaker variability. In this model GMM parameters are computed as a projection of the GMM supervector thus allowing to model linear combinations of speakers and each speaker has individual GMMs. 

Overall focus in state-of-art implementation has been shifted to DNNs for classification and they are better classifiers of course. Context-dependent tree to classify phonetic context is still constructed with GMMs though.