This is a brute force solution requiring just multivariable calculus.
It suffices to prove that the sample mean
$$
\bar{X} = \frac{1}{n} \sum_{i=1}^n X_i
$$
and the sample variance
$$
S^2 = \frac{1}{n - 1} \sum_{i=1}^n \left(X_i - \bar{X}\right)^2
$$
are independent. Thus, it suffices to prove that the sample mean $\bar{X}$ is independent of the vector
$$
(X_1 - \bar{X}, \ldots, X_n - \bar{X}).
$$
Moreover, since
$$
\begin{aligned}
\sum_{i=1}^n (X_i - \bar{X})
&= \sum_{i=1}^n X_i - \sum_{i=1}^n \bar{X} \\
&= n \bar{X} - n \bar{X} \\
&= 0,
\end{aligned}
$$
and hence
$$
X_1 - \bar{X}
= -\sum_{i=2}^n (X_2 - \bar{X}),
$$
it follows that $X_1 - \bar{X}$ can be recovered from just knowing $(X_2 - \bar{X}, \ldots, X_n - \bar{X})$.
Thus, it suffices to prove that the sample mean $\bar{X}$ is independent from
$$
(X_2 - \bar{X}, \ldots, X_n - \bar{X}).
$$
Now consider the joint density
$$
\begin{aligned}
f_{(X_1, \ldots, X_n)}(x_1, \ldots, x_n)
&= \left(2 \pi \sigma^2\right)^{-n/2} \exp\left(-\sum_{i=1}^n \frac{1}{2}\left(\frac{x_i - \mu}{\sigma}\right)^2\right) \\
&= \left(2 \pi \sigma^2\right)^{-n/2} \exp\left(-\sum_{i=1}^n \frac{1}{2}\left(\frac{x_i - \bar{x}}{\sigma}\right)^2 - \frac{n}{2}\left(\frac{\bar{x} - \mu}{\sigma}\right)^2\right) \\
&= \underbrace{\left(2 \pi \sigma^2\right)^{-n/2}}_{\text{constant}}
\underbrace{\exp\left(-\sum_{i=1}^n \frac{1}{2}\left(\frac{x_i - \bar{x}}{\sigma}\right)^2\right)}_{\text{depends only on $(x_2-\bar{x},\ldots,x_n-\bar{x})$}}
\underbrace{\exp\left(-\frac{n}{2}\left(\frac{\bar{x} - \mu}{\sigma}\right)^2\right)}_{\text{depends only on $\bar{x}$}}.
\end{aligned}
$$
To get from $(X_1,\ldots,X_n)$ to $(\bar{X}, X_2 - \bar{X}, \ldots, X_n - \bar{X})$, consider the diffeomorphism $T : \mathbb{R}^n \to \mathbb{R}^n$ given by
$$
T(x_1, \ldots, x_n)
= (\bar{x}, x_2 - \bar{x}, \ldots, x_n - \bar{x}).
$$
($T$ is a diffeomorphism since it's clearly differentiable and its inverse is given by
$$
T^{-1}(y_1, \ldots, y_n)
= \left(n y_1 - \sum_{i=2}^n y_i, y_2 + y_1, \ldots, y_n + y_1\right),
$$
which is also clearly differentiable).
Up to transpose, the Jacobian matrix of $T$ is
$$
DT(x_1, \ldots, x_n)
= \begin{bmatrix}
1/n & 1/n & 1/n & \cdots & 1/n \\
-1/n & (n - 1) / n & -1/n & \cdots & -1/n \\
-1/n & -1/n & (n - 1) / n & \cdots & -1/n \\
\vdots & \vdots & \vdots & \ddots & \vdots \\
-1/n & -1/n & -1/n & \cdots & (n - 1)/n.
\end{bmatrix},
$$
which doesn't depend on $x_1, \ldots, x_n$.
Thus, the determinant of $DT$ is some constant $C$.
Now the joint density of $(\bar{X}, X_2 - \bar{X}, \ldots, X_n - \bar{X})$ satisfies
$$
f_{(\bar{X}, X_2 - \bar{X}, \ldots, X_n - \bar{X})}(y_1, \ldots, y_n)
= |C| f_{(X_1, \ldots, X_n)}(T^{-1}(y_1, \ldots, y_n))
$$
which factors as a function of $y_1$ times a function of $(y_2, \ldots, y_n)$ by what was shown above.
Therefore, $\bar{X}$ and $(X_2 - \bar{X}, \ldots, X_n - \bar{X})$ are independent.
You can find a general result about the quadratic form of a normal random vector in this related question. In the case where the the normal random variables are independent wth common mean but different variances, the quadratic form will be a weighted sum of chi-squared-one random variables:
For the random variables specified in your question, you can group them into the normal random vector:
$$\mathbf{X} \sim \text{N}(\mu \mathbf{1}, \boldsymbol{\Sigma})
\quad \quad \quad \quad \quad
\boldsymbol{\Sigma} = \text{diag}(\sigma_1^2,...,\sigma_n^2).$$
Using the centering matrix $\mathbf{C}$ we have $\mathbf{X} - \bar{\mathbf{X}} = \mathbf{C} \mathbf{X}$. Moreover, it is simple to show that the matrix $\mathbf{C}$ is symmetric and idempotent, so that $\mathbf{C}^\text{T} \mathbf{C} = \mathbf{C}$. (For detailed information on the centering matrix, see e.g., O'Neill 2020, esp. sections 3-4.) For this analysis I'm also going to use the standard normal random vector $\mathbf{Z} \sim \text{N}(\mathbf{0}, \mathbf{I})$. Using these objects, we can write the quadratic form of interest as:
$$\begin{align}
\sum_{i=1}^n (X_i - \bar{X})^2
&= (\mathbf{X} - \bar{\mathbf{X}})^\text{T} (\mathbf{X} - \bar{\mathbf{X}}) \\[6pt]
&= (\mathbf{C} \mathbf{X})^\text{T} (\mathbf{C} \mathbf{X}) \\[6pt]
&= \mathbf{X}^\text{T} \mathbf{C}^\text{T} \mathbf{C} \mathbf{X} \\[6pt]
&= \mathbf{X}^\text{T} \mathbf{C} \mathbf{X} \\[6pt]
&= (\mathbf{X} - \mu \mathbf{1})^\text{T} \mathbf{C} (\mathbf{X} - \mu \mathbf{1}) \\[6pt]
&= (\boldsymbol{\Sigma}^{1/2} \mathbf{Z})^\text{T} \mathbf{C} (\boldsymbol{\Sigma}^{1/2} \mathbf{Z}) \\[6pt]
&= \mathbf{Z}^\text{T} \boldsymbol{\Sigma}^{1/2} \mathbf{C} \boldsymbol{\Sigma}^{1/2} \mathbf{Z} \\[6pt]
&\sim \sum_{i=1}^n \lambda_i \cdot \chi_1^2,
\end{align}$$
where $\lambda_1,...,\lambda_n$ are the eigenvalues of the matrix $\boldsymbol{\Sigma}^{1/2} \mathbf{C} \boldsymbol{\Sigma}^{1/2}$. This latter matrix is given by:
$$\begin{align}
\boldsymbol{\Sigma}^{1/2} \mathbf{C} \boldsymbol{\Sigma}^{1/2}
&= \text{diag}(\sigma_1,...\sigma_n) \ \mathbf{C} \ \text{diag}(\sigma_1,...\sigma_n) \\[6pt]
&= \begin{bmatrix}
\tfrac{n-1}{n} \sigma_1^2 & -\tfrac{1}{n} \sigma_1 \sigma_2 & \cdots & -\tfrac{1}{n} \sigma_1 \sigma_n \\
-\tfrac{1}{n} \sigma_1 \sigma_2 & \tfrac{n-1}{n} \sigma_2^2 & \cdots & -\tfrac{1}{n} \sigma_2 \sigma_n \\
\vdots & \vdots & \ddots & \vdots \\
-\tfrac{1}{n} \sigma_1 \sigma_n & -\tfrac{1}{n} \sigma_2 \sigma_n & \cdots & \tfrac{n-1}{n} \sigma_n^2 \\
\end{bmatrix}. \\[6pt]
\end{align}$$
You can take the above matrix form and use it to compute the eigenvalues $\lambda_1,...,\lambda_n$ for your standard deviation values $\sigma_1,...,\sigma_n$. This then gives you the weightings for the weighted sum of chi-squared-one random variables.
Best Answer
The distribution of $x_{n+1}$ conditional on $\overline x=(x_1 +\cdots +x_n)/n$ and $s^2 = \big((x_1-\overline x)^2+\cdots+(x_n-\overline x)^2\big)/(n-1)$ is $\operatorname N(\mu,\sigma^2),$ since the observations $x_1,\ldots,x_n,x_{n+1}$ are independent. Thus I suspect that what you actually want is the distribution of $(x_{n+1}-\overline x)/s,$ upon which you can base a prediction interval whose endpoints are $\overline x\pm c\cdot s_n,$ where $c$ is chosen so as to get a desired probability that $x_{n+1}$ is in the interval.
Notice that $x_{n+1}-\overline x\sim\operatorname N\left(0,\sigma^2\left( 1 + \frac 1 n \right)\right),$ and since $s^2$ is independent of $\overline x$ and of $x_{n+1},$ and $(n-1)s^2/\sigma^2\sim\chi^2_{n-1},$ you have $$ \frac{(x_{n+1}-\overline x)/\sqrt{1+\frac 1 n}}{s/\sqrt n} \sim t_{n-1}. $$ Thus if you choose $c$ so that $\Pr(-c<t_{n-1}<+c)= 1-\alpha,$ then you have $$ \Pr\left(x_{n+1} \text{ is between } \overline x \pm c\cdot\frac s {\sqrt n}\cdot\sqrt{1+\tfrac 1 n} \right) = 1-\alpha. $$ This is a prediction interval for $x_{n+1}.$