This question resembles Vasile Cirtoaje’s equal-variables results, as explained here (although those results may be useless in the present problem).
Let $s,p$ be two positive numbers with $\frac{s}{3} \gt p^{\frac{1}{3}}$, and let
$$ K=\lbrace (x,y,z) \in {\mathbb R}^3 \big| 0 \lt x \leq y \leq z, x+y+z=s, xyz=p \rbrace $$
The set $K$ is compact, so the function $f(x,y,z)=(1+x^2)(1+y^2)(1+z^2)$ attains its maximum on $K$ at some point $(x_1,y_1,z_1) \in K$. Is it always true that $x_1=y_1$ or $y_1=z_1$ ?
Note that both cases can occur. For example, if $s=6$ and $p=7$ then $y_1=z_1$ but if $s=\frac{343}{180}$ and $p=\frac{171}{2000}$ then $x_1=y_1$.
It is straightforward enough to see that $K$ can be written
$$ K=\lbrace (x,\frac{s-x-\sqrt{(s-x)^2-\frac{4p}{x}}}{2},\frac{s-x+\sqrt{(s-x)^2-\frac{4p}{x}}}{2}) | x \in [\alpha,\beta]\rbrace $$
where $\alpha$ is the smallest root of $A(x)=x(s-x)^2-4p$, and $\beta$ is the second smallest root of $B(x)=x(s-3x)^2-A(x)$. The question can then be rephrased as, is it true that
$$ g(x)=f (x,\frac{s-x-\sqrt{(s-x)^2-\frac{4p}{x}}}{2},\frac{s-x+\sqrt{(s-x)^2-\frac{4p}{x}}}{2}) $$
always attains its maximum on the boundary of its domain of definition, at $x=\alpha$ or $x=\beta$.
Putative counterexamples may be found by looking at what happens when $g(\alpha)=g(\beta)$. When $s=\frac{343}{180}$, for example, there is a unique value $p_0$ of $p$ (with $p \approx 0.1385 \ldots$) making $g(\alpha)=g(\beta)$. It is obvious that $p_0$ is an algebraic number, and I’ve computed $p_0$ to 100 decimal places, but that did not help in finding the degree of $p_0$ (which seems to be prohibitively large).
It can also be shown that if $(x_1,y_1,z_1)\in K$ is a local maximum for $f$ (restricted to $K$) and $x_1 \neq y_1, y_1 \neq z_1$, then $x_1y_1+x_1z_1+y_1z_1=1$ (this follows from using the implicit function theorem and a local expansion around $(x_1,y_1,z_1)$ for $f$).
It is true.
The key idea is to reparametrize by $x+y+z$, $xy+yz+zx$, $xyz$ instead, and note that your function $$f(x,y,z) = (xy+yz+zx)^2 - 2(xy+yz+zx) + constant$$ when we fix $x+y+z$ and $xyz$. The problem is then about the restriction of range under this change of parametrization. Since $f$ is convex in terms of $xy+yz+zx$, it suffices to have a result that says $xy+yz+zx$ reaches its maximum/minimum with $x+y+z$,$xyz$ fixed when two of the variables are equal. This is already studied, under the name of uvw method. See here for example, "Tej's theorem" part II on p.3.
I guess this is what you are getting at. This does seem to reduce that problem to a single variable one.