The full answer to your question of why the answer is still normal requires appeal to a highly desirable property of the normal distribution: that it is a stable distribution. The theory of stable distributions requires characteristic equations and somewhat advanced, measure theoretic/real analytic arguments. Suffice it to say, most statisticians simply accept this fact, and henceforth just remember that linear combinations of normal variables are indeed normal.
In your case, your formulation is essentialy creating a multivariate normal vector with a gaussian "effect" added to it. Therefore, the resulting multivariate vector will again be a multivariate normal, but with a mean vector of $\mathbf{1}$ (since each X has mean 0 and Z has mean 1). However, the presence of the same Z in each term means that the resulting multivariate distribution is no longer a set of iid variables, since the common Z induces a correlation. As Dilip has pointed out, these covariances will be:
$cov(X_i+Z,X_j+Z)=cov(X_i,X_j)+cov(X_i,Z)+cov(Z,X_j)+cov(Z,Z)=\sigma^2_Z$ As all the other terms become zero due to independence.
Therefore, the correlation between any two elements of the resulting multivariate normal vector, ($X_i+Z,X_j+Z)$, will be $\rho = \frac{cov(X_i+Z,X_j+Z)}{\sqrt{\sigma_i^2+\sigma_Z^2}+ \sqrt{\sigma_j^2+\sigma_Z^2}}=\frac{4}{2\sqrt{5}}\approx .0.89$ Resulting in quite a strong linear relationship between the variables.
$\displaystyle \sum_{i=1}^N \left(\frac{x_i-\mu}{\sigma}\right)^2$ has a $\chi_N^2$ distribution (chi-squared with $N$ degrees of freedom) as as the sum of $N$ independent standard normal random variables.
So if $\mu=0$ and $\displaystyle s_2^2=\frac{1}{N}\sum_{i=1}^Nx_i^2$ then $N s_2^2 / \sigma^2$ also has a $\chi_N^2$ distribution.
Best Answer
The random variable $Z = \max_{i=1}^n(X_i)$ is known as order statistics, and is sometimes denoted as $X_{n:n}$.
The cumulative density function of $Z$ is easy to find: $$ F_Z(z) = \mathbb{P}\left(Z \leqslant z\right) = \mathbb{P}\left(\max_{i=1}^n(X_i) \leqslant z\right) = \mathbb{P}\left( X_1 \leqslant z, X_2 \leqslant z, \ldots, X_n \leqslant z\right) $$ using independence: $$ F_Z(z) = \left(F_X(z)\right)^n $$ Thus the density function is $$ f_Z(z) = n f_X(z) F_X^{n-1}(z) $$ In particular, it follows that $Z$ is not normal.
Expected values of $Z$ are known in closed form for $n=1,2,3,4,5$ (asking Mathematica):
Large $n$ asymptotics is discussed here.