Are linear combinations (i.e. "sums") of gaussian distributions also gaussian? I have always assumed that this is fact : but are there any mathematical theorems that state: finite sums of gaussian distributions are gaussian themselves?

In the above picture, does p(x) have a gaussian distribution? Do the sums need to be linear? Can they be non-linear?
I saw the following question: Is the joint distribution of two linear combinations of Gaussians still a multivariate normal? - but no reference/proof was provided.
Thanks

It is the linear combination of jointly Gaussian random variables (RVs) that results in another RV with Gaussian density. In your question, you have linear combination of Gaussian densities; therefore, the resulting density need not be Gaussian.
Below is given a working proof of this theorem. The characteristic function of an RV $X$ is $$\phi_X(t)=E(e^{\iota tx}),$$ where $E$ denotes expectation. The characteristic function of the linear combination $Y$ of two RVs $X_1$ and $X_2$ is $$Y=\alpha_1X_1+\alpha_2X_2,$$ and the characteristic function of $Y$ is $$\phi_Y(t)=E(e^{\iota tY})=E(e^{\iota t\alpha_1X_1 +\iota t\alpha_2X_2}).$$ If $X_1$ and $X_2$ are independent $\phi_Y(t)$ can be written as $$\phi_Y(t) = E(e^{\iota t\alpha_1X_1})E(e^{\iota t\alpha_2X_2}).$$
Now, for a Gaussian RV $X$ with mean $\mu$ and variance $\sigma^2$,$$\phi_X(t) = \exp(\iota t\mu-\frac{1}{2}\sigma^2t^2).$$ It is a bit tideous to derive this result but it can be done with some effort. Following this result, $$\phi_Y(t)=\exp[\iota t(\alpha_1\mu_1+\alpha_2\mu_2)-\frac{1}{2}(\alpha_1^2\sigma_1^2+\alpha_2^2\sigma_2^2)t^2].$$
Now, we also know that characteristic function of a probability distribution is unique (but it is difficult to prove). Following this result, characteristic function of $Y$ is clearly that of a Gaussian with mean $\alpha_1\mu_1 + \alpha_2\mu_2$ and variance $\alpha_1^2\sigma_1^2 + \alpha_2^2\sigma_2^2$. It basically proves that linear combination of two independent Gaussian RVs is a Gaussian RV.
Edit: The proof above is valid only when $X_1$ and $X_2$ are independent. But a similar proof can be carried out when $X_1$ and $X_2$ are not independent. Suppose $X = [X_1,X_2]^T$ is a jointly Gaussian with mean vector $\mu = [\mu_1,\mu_2]^T $ and covariance matrix $C$. Also, the characteristic function of a random vector $X$ is $$\phi_X(t_1,t_2)=E(e^{\iota t^TX})=E(e^{\iota (t_1X_1+t_2X_2)}),$$
where $t=[t_1,t_2]^T$. The characteristic function of an RV $Y=\alpha_1X_1+\alpha_2X_2$ is $$\phi_Y(u)=E(e^{\iota uY})=E(e^{\iota u(\alpha_1X_1+\alpha_2X_2)})$$ $$=E(e^{\iota (u\alpha_1X_1+u\alpha_2X_2)}).$$
Notice that, $$\phi_Y(u)=\phi_X(\alpha_1u,\alpha_2u).$$
Now, we know that (the proof is a bit involved): $$\phi_X(t)=\exp[\iota t^T\mu-\frac{1}{2}t^TCt],$$ which implies
$$\phi_Y(u)=\exp[\iota u\mu_y-\frac{1}{2}u^2\sigma_y^2],$$
where $$\mu_y=\alpha^T\mu=\alpha_1\mu_1+\alpha_2\mu_2,$$ $$\sigma_y^2=\alpha^TC\alpha$$ $$\alpha=[\alpha_1,\alpha_2]^T.$$
The last four equations show that $Y$ is a Gaussian RV with mean $\mu_y$ and variance $\sigma_y^2$. I have not given the derivation of the last four equations, let me know if that would be helpful.
A side note: The RVs $X_1$ and $X_2$ should be jointly Gaussian, otherwise this result would not hold. An illustrative example is given by Chris Huang in comments section.
