I suspect the only symbolic results you'll be able to find involve those distributions whose sum of iid random variables follow the same family of distributions. That is certainly the case for normal and gamma distributions.
Rather than show a step-by-step calculus approach, I'm going to just give handwaving arguments and Mathematica commands and I've slightly differed from your notation. You can fill in details if desired.
Normal distributions
If $X_i \sim N(\mu,\sigma^2)$, then $Y=\sum_{i=2}^n X_k$ has a normal distribution with mean $(n-1)\mu$ and variance $(n-1)\sigma^2$. $X_1$ and $S_n=X_1+Y$ have a bivariate normal distribution with mean vector $(\mu, n\mu)$ and covariance matrix
$$\Sigma=\left(
\begin{array}{cc}
\sigma ^2 & \sigma ^2 \\
\sigma ^2 & n \sigma ^2 \\
\end{array}
\right)$$
The conditional density of $X_1$ given $S_n=s$ is found by dividing the joint density of $(X_1, S_n))$ by the marginal density of $S_n$. One then integrates $x_1^2$ times the conditional density over values of $x_1$.
Joint distribution of $X_1$ and $S_n$.
dist = TransformedDistribution[{x1, x1 + x2n},
{x1 \[Distributed] NormalDistribution[μ, σ],
x2n \[Distributed] NormalDistribution[(n - 1) μ, σ Sqrt[n - 1]]}];
(* PDF's of joint and marginal distributions *)
pdfx1S = PDF[dist, {x1, s}]
$$\frac{\exp \left(\frac{1}{2} \left(-\frac{(\text{x1}-\mu ) (n \text{x1}-s)}{(n-1) \sigma ^2}-\frac{(-\mu -\mu (n-1)+s) (\mu -\mu n+s-\text{x1})}{(n-1) \sigma ^2}\right)\right)}{2 \pi \sqrt{n \sigma ^4-\sigma ^4}}$$
pdfs = Simplify[PDF[MarginalDistribution[dist, 2], s], Assumptions -> s > 0]
$$\frac{e^{-\frac{(s-\mu n)^2}{2 n \sigma ^2}}}{\sqrt{2 \pi } \sqrt{n \sigma ^2}}$$
(* Find conditional expectation by integrating *)
Integrate[x1^2 pdfx1S/pdfs, {x1, -∞, ∞}, Assumptions -> {σ > 0, n > 1}]
$$\frac{(n-1) n \sigma ^2+s^2}{n^2}$$
Gamma distributions
If $X_i \sim \Gamma(k, \theta)$, then $Y=\sum_{i=2}^n X_i$ has a gamma distribution with parameters $(n-1)k$ and $\theta$. We follow the same approach as above.
Joint distribution of $X_1$ and $S_n=X_1+Y$.
dist = TransformedDistribution[{x1, x1 + x2n},
{x1 \[Distributed] GammaDistribution[k, θ],
x2n \[Distributed] GammaDistribution[(n - 1) k, θ]}];
(* PDF's of joint and marginal distributions *)
pdfx1S = Simplify[PDF[dist, {x1, s}] // PiecewiseExpand, Assumptions -> {s > x1 > 0}]
$$\frac{\text{x1}^{k-1} e^{-\frac{s}{\theta }} \theta ^{-k n} (s-\text{x1})^{k (n-1)-1}}{\Gamma (k) \Gamma (k (n-1))}$$
pdfs = Simplify[PDF[MarginalDistribution[dist, 2], s], Assumptions -> s > 0]
$$\frac{e^{-\frac{s}{\theta }} \theta ^{-k n} s^{k n-1}}{\Gamma (k n)}$$
(* Find conditional expectation by integrating *)
Integrate[x1^2 pdfx1S/pdfs, {x1, 0, s}, Assumptions -> {s > 0, k > 0, n > 1}]
$$\frac{(k+1) s^2}{n (k n+1)}$$
Conclusion
In short, the condition expectation depends on the common distribution and it is unlikely to obtain a symbolic solution for very many distribution families. It might seem odd that for the normal the conditional expectation is not a function of $\mu$ and for the gamma distribution the conditional expectation is not a function of $\theta$.
