Take $t = 5k$ and $z = 7k$, then we want the equality $x^2 + y^2 = k^2$ to hold. Consider two series of solutions: $x = n^2 - 1$, $y = 2n$, $k = n^2 + 1$ and $x = 3n^2 + 4n + 1$, $y = 2n(2n+1)$, $k = 5n^2 + 4n + 1$. It is evident that in both cases $(x, y, z, t) = 1$ (here and on $(x, y, z, t)$ denotes $\text{GCD}\{x, y, z, t\}$) for even $n$. Consider the first series and denote $Q = 2P^2(t) + 2P(xy + yz + zx) - P^2(x + y + z)$ and $P(x) = a_m x^m + a_{m-1} x^{m-1} + \dots + a_1 x + a_0$, $a_m \neq 0$. Then polynomial $Q_1(n) = 2P^2(5(n^2+1)) + 2P(7n^4 + 16n^3 - 12n - 7) - P^2(8n^2 + 2n + 6)$ equals to zero for every $n = 2s$ with $s \in \mathbb{N}$, whence $Q_1(n) = 0$. Consider the coefficient at $x^{4m}$ in $Q_1(x)$: on the one hand it is equal to $2 \cdot 5^{2m} a_m^2 + 2 \cdot 7^m a_m - 8^{2m} a_m^2$ and on the other $0$, therefore $2 \cdot 5^{2m} a_m^2 + 2 \cdot 7^m a_m - 8^{2m} a_m^2 = 0$. Analogously we can obtain that $Q_2(n) = 2P^2(5(5n^2+4n+1)) + 2P(257n^4 + T_2(n)) - P^2(42n^2 + 34n + 8)$, with $\deg(T_2) \le 3$, is equal to zero and thus the coefficient at $x^{4m}$ in this polynomial also equals zero: $2 \cdot 25^{2m} a_m^2 + 2 \cdot 257^m a_m - 42^{2m} a_m^2 = 0$. So, $a_m = \frac{2 \cdot 7^m}{64^m - 2 \cdot 25^m} = \frac{2 \cdot 257^m}{42^{2m} - 2 \cdot 25^{2m}}$, whence $2 \cdot 7^m \cdot (42^{2m} - 2 \cdot 25^{2m}) = (64^m - 2 \cdot 25^m) \cdot 2 \cdot 257^m$. For $m \ge 2$ we have: $(64^m - 2 \cdot 25^m) \cdot 2 \cdot 257^m > (64^m - 2 \cdot 25^m) \cdot 2 \cdot 252^m$, and so $42^{2m} - 2 \cdot 25^{2m} > 36^m (64^m - 2 \cdot 25^{2m})$. $7^{2m} > 7^{2m} - \frac{2 \cdot 25^{2m}}{36^m} > 64^m - 2 \cdot 25^m = (7^2 + 15)^m - 2 \cdot 25^m > 7^{2m} + m \cdot 7^{2m-2} \cdot 15 - 2 \cdot 25^m \Rightarrow 2 \cdot 25^m > m \cdot 49^{m-1} \cdot 15 \ge 2 \cdot 49^{m-1} \cdot 15 \Rightarrow \frac{49^{m-1}}{15} \ge \frac{49}{25}$, which is evidently false. Thus, $m \le 1$. For $m = 1$, we have $P(x) = a_1 x + a_0$ and $a_1 = 1$, i.e. $P(x) = x + a_0$. Substitute $P(x)$ into the initial equality and take $(x, y, z, t)$ from the first series: $(n^2 - 1, 2n, 7(n^2 + 1), 5(n^2 + 1))$. After some easy transformations we will have: $4a_0(5n^2 + 5) + 2a_0^2 + 2a_0 - 2a_0(8n^2 + 2n + 6) - a_0^2 = 0$, which implies that $a_0 = 0$ or that $a_0$ depends on $n$, which is impossible. So for $m = 1$ we have solution $P(x) = x$. Now if $n = 0$ then $a_0^2 + 2a_0 = 0$, and we obtain two more answers: $P(x) = 0$ and $P(x) = -2$. An easy check shows that these solutions meet all the requirements.

Consider the equation $x^2 + y^2 + z^2 = 2t^2$ and make use of the method of Diophantine secants: $\left(\frac{x}{t}\right)^2 + \left(\frac{y}{t}\right)^2 + \left(\frac{z}{t}\right)^2 = 2$, i.e. $\alpha^2 + \beta^2 + \gamma^2 = 2$, where $\alpha, \beta, \gamma \in \mathbb{Q}^+$. Evidently, $\alpha = 0$, $\beta = \gamma = 1$ is a solution of the last equation. Now take $\beta = -k_1\alpha + 1$, $\gamma = -k_2\alpha + 1$, with $k_1, k_2 \in \mathbb{Q}$. Then $\alpha^2 + (1-k_1\alpha)^2 + (1-k_2\alpha)^2 = 2 \Leftrightarrow \alpha(k_1^2\alpha + k_2^2\alpha - 2k_1 - 2k_2 + \alpha) = 0$. If $\alpha \neq 0$, we have $\alpha = \frac{2(k_1+k_2)}{1+k_1^2+k_2^2} \Rightarrow \beta = \frac{1+k_2^2-k_1^2-2k_1k_2}{1+k_1^2+k_2^2}$, $\gamma = \frac{1+k_1^2-k_2^2-2k_1k_2}{1+k_1^2+k_2^2}$.
Take $k_1 = \frac{n}{p}$, $k_2 = \frac{m}{p}$, $n, m, p \in \mathbb{N}$. Then $\alpha = \frac{2(m+n)p}{p^2+n^2+m^2}$, $\beta = \frac{p^2+m^2-n^2-2mn}{p^2+n^2+m^2}$, $\gamma = \frac{p^2+n^2-m^2-2mn}{p^2+n^2+m^2}$. Thus, we see that $x = 2(m+n)p$, $y = p^2+m^2-n^2-2mn$, $z = p^2+n^2-m^2-2mn$, $t = p^2+m^2+n^2$, with $n, m, p \in \mathbb{N}$ are the solutions of the equation $x^2 + y^2 + z^2 = 2t^2$.
First consider $m = 2$, $n = 1$. Then $x_1 = 6p$, $y_1 = p^2-1$, $z_1 = p^2-7$, $t_1 = p^2+5$. It is evident that for $p \neq 6$ we have $(x_1, y_1, z_1, t_1) = 1$. Then substituting these $(x_1, y_1, z_1, t_1)$ into the problem statement we obtain that $Q = 2P^2(t) + 2P(xy + yz + zx) - P^2(x + y + z) = 0$, which means that $Q_1(p) = 2P^2(p^2+5) + 2P(p^4+12p^3-8p^2-48p+7) - P^2(2p^2+6p-8) = 0$. From the problem conditions it follows that all numbers $p$ of the form $p=6k$ should be the roots of $Q_1$, thus $Q_1 \equiv 0$.
Let $P(x) = a_n x^n + a_{n-1} x^{n-1} + \dots + a_1 x + a_0$, with $a_n \neq 0$. Then $Q_1(p) = 2(a_n(p^2+5)^n + \dots + a_0)^2 + 2(a_n(p^4+12p^3-8p^2-48p+7)^n + \dots + a_0)^2 - (a_n(2p^2+6p-8)^n + \dots + a_0) = (2a_n^2 + 2a_n - 2^{2n}a_n^2)p^{4n} + R_1(p)$, where $\deg R_1 \le 4n-1$. Since $Q_1(p) \equiv 0$ and $a_n \neq 0$ we have that
$$2a_n + 2 - 2^{2n}a_n = 0. \quad (1)$$
Consider $p = 2m$, $n = 1$. Then $x_2 = 4m(m+1)$, $y_2 = 5m^2 - 2m - 1$, $z_2 = 3m^2 - 2m + 1$, $t_2 = 5m^2 + 1$. For $m$ large enough we have $t_2 > y_2 > x_2 > z_2$, and so these numbers are distinct. If $m$ is even, $(x_2, z_2) = (4m(m+1), 3m^2 - 2m + 1) = (4(m+1), 3m^2 - 2m + 1) = (m+1, 3m^2 - 2m + 1) = (m+1, 5m-1) = (m+1, 6)$, thus for $m \neq 6$ $(x_2, y_2, z_2, t_2) = 1$. If we follow the same lines as for the first series $(x, y, z, t)$, we will obtain that $Q_2(m) = 2P^2(5m^2+1) + 2P(47m^4 + T(m)) - P^2(12m^2) = 0$, where $\deg T \le 3$. Let us compute the coefficient at $m^{4n}$ in $Q_2(m)$ - it's equal to $2a_n^2 \cdot 5^{2n} + 2a_n \cdot 47^n - 12^{2n}a_n^2 = 0$. Since $a_n \ne 0$ we have that:
$$2a_n \cdot 5^{2n} + 2 \cdot 47^n - 12^{2n} a_n = 0. \qquad (2)$$
Now from (1) and (2) we obtain that $a_n = \frac{1}{2^{2n-1}}$ and $(2^{2n-1}-1) \cdot 2 \cdot 47^n = 12^{2n} - 2 \cdot 5^{2n}$, whence $2 \cdot 5^{2n} + (4 \cdot 47)^n = 144^n + 2 \cdot 47^n$. It's easy to see that the last equality is impossible for arbitrary natural $n \ge 2$, since $2 \cdot 5^{2n} + (4 \cdot 47)^n > 188^n = (144 + 44)^n \ge 144^n + 44 \cdot 144^{n-1} > 144^n + 2 \cdot 47^n$.
So $\deg P \le 1$. First assume that $\deg P = 1$. Then the equality $a_n = \frac{1}{2^{2n-1}}$ implies that $a_1 = 1$, i.e. that $P(x) = x + a_0$. Substitute this polynomial into the initial equality:
$$2P^2(t) + 2P(xy + yz + zx) = P^2(x + y + z) \Rightarrow \\ 2(t^2 + 2a_0t + a_0^2) + 2(xy + yz + zx) + 2a_0 = (x + y + z)^2 + 2(x + y + z)a_0 + a_0^2 \Leftrightarrow \\ 4a_0t + a_0^2 + 2a_0 = 2(x + y + z)a_0.$$
If $a_0 = 0$ then $P(x) = x$, otherwise $4t + a_0 + 2 = 2(x + y + z)$. Substitute here the known series $x_1 = 6p$, $y_1 = p^2 - 1$, $z_1 = p^2 - 7$, $t_1 = p^2 + 5$ with $p = 6$. Then $4(p^2 + 5) + a_0 + 2 = 2(2p^2 + 6p - 8)$, and we get that the constant $a_0$ depends on $p$, which is impossible.
If $\deg P = 0$ then $P = a_0 = \text{const}$. From the problem conditions we obtain that $2a_0^2 + 2a_0 = a_0^2 \Rightarrow a_0 = -2$ also $a_0 = 0$, which means that $P(x) = 0$ or $P(x) = -2$. An easy check shows that all three polynomials obtained meet all the requirements.