How to show that any Gaussian time-series is linear one? In this paper I saw the following statement:

If the time series is Gaussian (i.e., normally distributed) then the
  best linear forecast is in fact the best of all possible forecasts: No
  nonlinear forecast can do better in terms of mean squared prediction
  error. Thus, as long as the series is Gaussian, we need look no
  further than the linear methods (e.g., ARMA forecasting) already
  presented.

I think it is a very strong, general and interesting statement. Is it easy to prove, that a time-series whose values are normally distributed are given ("generated") by a linear model? Or, alternatively, any nonlinear time-series should have a non Gaussian distribution of values.
 A: This is more general than the normal case.  
If we take as criterion of optimality the minimization of Mean Squared Forecasting Error (and it is not the only one, but perhaps it is the more widely used), then the best predictor is the Conditional Expectation (i.e. the conditional expected value treated as a function), for any distributional assumption. There are some threads here in CV that contain this proof, for example, https://stats.stackexchange.com/a/71869/28746
Then, for those joint distributions  that result in Conditional Expected Value being a linear function of the conditioning variables, it follows that the best in minimum MSFE sense predictor will be a linear function of the data.
Linear Conditional Expectations characterize the symmetric elliptical class as well as many of the Pearson family of multivariate distributions. The normal distribution belongs here.
A: 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.
