[proofplan]
We identify the distribution of the sum by computing its characteristic function. The characteristic function of each Gaussian summand has the exponential form determined by its mean and variance, and independence lets us multiply these characteristic functions. The resulting product is exactly the characteristic function of a [Gaussian random variable](/page/Gaussian%20Random%20Variable) whose mean and variance are the corresponding sums. The uniqueness theorem for characteristic functions then gives equality in distribution.
[/proofplan]
[step:Introduce the sum and its characteristic function]
Define the [random variable](/page/Random%20Variable) $S:(\Omega,\mathcal F)\to(\mathbb R,\mathcal B(\mathbb R))$ by
\begin{align*}
S(\omega)=\sum_{j=1}^n X_j(\omega).
\end{align*}
For each real-valued random variable $Y:(\Omega,\mathcal F)\to(\mathbb R,\mathcal B(\mathbb R))$, define its characteristic function $\phi_Y:\mathbb R\to\mathbb C$ by
\begin{align*}
\phi_Y(t)=\mathbb E[e^{itY}].
\end{align*}
The random variable $e^{itY}$ is bounded by $1$ in absolute value, so this expectation is finite for every $t\in\mathbb R$.
[/step]
[step:Compute the characteristic function of the sum by independence]
Fix $t\in\mathbb R$. Since $X_1,\dots,X_n$ are independent, the random variables $e^{itX_1},\dots,e^{itX_n}$ are independent and bounded. Hence the standard factorization formula for expectations of products of bounded independent random variables gives
\begin{align*}
\phi_S(t)=\mathbb E[e^{itS}]=\mathbb E\left[\prod_{j=1}^n e^{itX_j}\right]=\prod_{j=1}^n \mathbb E[e^{itX_j}]=\prod_{j=1}^n \phi_{X_j}(t).
\end{align*}
Using the standard formula for the characteristic function of a Gaussian random variable, for each $j\in\{1,\dots,n\}$ and each $t\in\mathbb R$,
\begin{align*}
\phi_{X_j}(t)=\exp\left(i\mu_j t-\frac{\sigma_j^2t^2}{2}\right).
\end{align*}
Therefore
\begin{align*}
\phi_S(t)=\prod_{j=1}^n \exp\left(i\mu_j t-\frac{\sigma_j^2t^2}{2}\right).
\end{align*}
Combining the exponents gives
\begin{align*}
\phi_S(t)=\exp\left(it\sum_{j=1}^n\mu_j-\frac{t^2}{2}\sum_{j=1}^n\sigma_j^2\right).
\end{align*}
[guided]
The goal is to convert the independence of the random variables into an algebraic identity for the characteristic function of their sum. Fix $t\in\mathbb R$. By definition,
\begin{align*}
\phi_S(t)=\mathbb E[e^{itS}].
\end{align*}
Since $S=\sum_{j=1}^n X_j$, the exponential law gives
\begin{align*}
e^{itS}=\prod_{j=1}^n e^{itX_j}.
\end{align*}
The maps $x\mapsto e^{itx}$ are Borel measurable from $\mathbb R$ to $\mathbb C$, so independence of $X_1,\dots,X_n$ implies independence of $e^{itX_1},\dots,e^{itX_n}$. Each factor satisfies $|e^{itX_j}|=1$, so the factors are bounded and integrable. Thus the factorization formula for bounded independent random variables applies:
\begin{align*}
\phi_S(t)=\mathbb E\left[\prod_{j=1}^n e^{itX_j}\right]=\prod_{j=1}^n\mathbb E[e^{itX_j}]=\prod_{j=1}^n\phi_{X_j}(t).
\end{align*}
Now we insert the Gaussian characteristic function formula. Since $X_j\sim\mathcal N(\mu_j,\sigma_j^2)$, the characteristic function of $X_j$ is
\begin{align*}
\phi_{X_j}(t)=\exp\left(i\mu_j t-\frac{\sigma_j^2t^2}{2}\right).
\end{align*}
Therefore
\begin{align*}
\phi_S(t)=\prod_{j=1}^n \exp\left(i\mu_j t-\frac{\sigma_j^2t^2}{2}\right).
\end{align*}
The product of exponentials is the exponential of the sum of the exponents, so
\begin{align*}
\phi_S(t)=\exp\left(\sum_{j=1}^n i\mu_jt-\sum_{j=1}^n\frac{\sigma_j^2t^2}{2}\right).
\end{align*}
Since $t$ is fixed and finite sums distribute over scalar multiplication, this becomes
\begin{align*}
\phi_S(t)=\exp\left(it\sum_{j=1}^n\mu_j-\frac{t^2}{2}\sum_{j=1}^n\sigma_j^2\right).
\end{align*}
This is the exact place where independence is used: without it, the characteristic function of the sum would not generally factor into the product of the characteristic functions.
[/guided]
[/step]
[step:Identify the computed characteristic function as Gaussian]
Define
\begin{align*}
m=\sum_{j=1}^n\mu_j
\end{align*}
and
\begin{align*}
v=\sum_{j=1}^n\sigma_j^2.
\end{align*}
Since each $\sigma_j^2\geq 0$, we have $v\geq 0$. Let $G$ be a Gaussian random variable with distribution $\mathcal N(m,v)$. By the Gaussian characteristic function formula, its characteristic function $\phi_G:\mathbb R\to\mathbb C$ satisfies
\begin{align*}
\phi_G(t)=\exp\left(im t-\frac{v t^2}{2}\right).
\end{align*}
Substituting the definitions of $m$ and $v$ gives, for every $t\in\mathbb R$,
\begin{align*}
\phi_G(t)=\exp\left(it\sum_{j=1}^n\mu_j-\frac{t^2}{2}\sum_{j=1}^n\sigma_j^2\right)=\phi_S(t).
\end{align*}
[/step]
[step:Apply uniqueness of characteristic functions]
The characteristic functions $\phi_S$ and $\phi_G$ agree at every $t\in\mathbb R$. By the uniqueness theorem for characteristic functions, equality of characteristic functions on $\mathbb R$ implies equality in distribution. Hence $S$ and $G$ have the same distribution, so
\begin{align*}
\sum_{j=1}^n X_j=S\sim\mathcal N\left(\sum_{j=1}^n\mu_j,\sum_{j=1}^n\sigma_j^2\right).
\end{align*}
This is the desired conclusion.
[/step]