In my opinion, there are two different questions being asked here, and the answer to the first is a Yes.
Given that a time series is Gaussian (meaning that all choices
of $n \geq 1$ random variables $X_{t_1}, X_{t_2}, \cdots, X_{t_n}$ from
the processes are jointly Gaussian random variables), is it true
that the minimum mean-square-error (MMSE) predictor of $X_N$ in terms of
any other variables from the series is a linear (or affine) function
of the known variables?
As Alecos's answer points out, this is a property enjoyed by
jointly Gaussian random variables. It is a very general
result (applicable to all random variables, not just Gaussian or
jointly Gaussian ones) that the MMSE predictor
is the conditional mean of $X_N$ given the other random variables.
A simpler predictor of $X_N$ is a linear combination (technically
an affine function) of the known variables, and if this function
has the smallest mean-square error among all possible linear
predictors, then it is called the Linear MMSE (LMMSE) predictor
of $X_N$.
For jointly Gaussian random variables, the conditional mean
(i.e. MMSE predictor) is
a linear function of the known variables, and so is also the
LMMSE predictor. Thus, the methods (already known to the OP)
that produce LMMSE estimators cannot be
improved upon for Gaussian time series: those LMMSE estimators
are in fact the MMSE estimators, and so we need not have any
niggling doubts that had we put in more hard work, we might
have been able to come up with a nonlinear predictor that is
better than the "easily-found" LMMSE estimator even if the
said nonlinear predictor is not as good as the MMSE estimator.
It is worthwhile noting that while most other random variables
do not enjoy the property that the LLMSE estimator and
the MMSE estimator coincide, it is by no means true that
jointly Gaussian random variables are unique in this respect.
Consider, for example $(X,Y)$ being uniformly distributed
on interior of the triangle with vertices $(0,0), (1,1), (1,0)$.
The conditional density of $Y$ given $X=x$ is $U(0,x)$ and
thus has mean $\frac x2$. Thus, the MMSE estimator of
$Y$ given $X$ is $\frac X2$ which is a linear function of
$X$ and hence must be the same as the LMMSE estimator of
$Y$ given $X$.
The other question being asked seems to be
Are all Gaussian time series generated by a "linear model"?
or
Are all time series generated by a "linear model" Gaussian time series?
or
Gaussian series cannot be generated by nonlinear models
If by linear model is meant that each random variable $X_N$
in the series can be
expressed as
$$X_N = a_N + \sum_{i} b_{i,N} Y_i$$
where the $a$'s and $b$'s are constants
and $\{Y_i\}$ is another time series, then the answer to the
middle question is No unless of course $\{Y_i\}$ itself is a
Gaussian time series. Indeed, since any collection of $L$
jointly Gaussian random variables is obtained via a linear
transformation of $M \leq L$ independent standard normal random
variables, the answer to the first question is Yes. On
the third question, I maintain a discreet silence.
