This is not an answer, but a reduction to supposedly simpler problem. (edited with more details)

Using quadratic Gauss sums, we can express Legendre symbol $\left(\frac{i^2+j^2}p\right)$ as \begin{split} \left(\frac{i^2+j^2}p\right) &= \frac{1}{I\sqrt{p}}\sum_{n=0}^{p-1} \zeta^{(i^2+j^2)n^2} \\ &= \frac{1}{I\sqrt{p}}\left(1+2\sum_{n=1}^{(p-1)/2} \zeta^{i^2n^2}\zeta^{j^2n^2}\right),\qquad(\star) \end{split} where $\zeta := \exp(\frac{2\pi I}p)$ and $I$ is the imaginary unit. From $(\star)$ it follows that $$\left[\left(\frac{i^2+j^2}p\right)\right]_{i,j=1}^{(p-1)/2} = \frac{1}{I\sqrt{p}}(J + 2B^2),$$ where $J$ is the matrix of all ones and $B:=[\zeta^{i^2j^2}]_{i,j=1}^{(p-1)/2}$. The enables application of pavl0's approach in evaluating the determinant in order to show that it is a square.

Alternatively, it is possible to "extract" the square root from the determinant without knowing its value. To do so, we need to embed "$1+$" into the quadratic sum in $(\star)$ by introducing a parameter $\alpha$: $$1+2\sum_{n=1}^{(p-1)/2} \zeta^{i^2n^2}\zeta^{j^2n^2} = 2\sum_{n=1}^{(p-1)/2} (\zeta^{i^2n^2}+\alpha)(\zeta^{j^2n^2}+\alpha).$$ Since $2\sum_{n=1}^{(p-1)/2} \zeta^{i^2n^2} = I\sqrt{p} - 1$ for any $i\not\equiv 0\pmod{p}$, we get that $\alpha$ satisfies the quadratic equation: $$1 = 2(I\sqrt{p} - 1)\alpha + (p-1)\alpha^2.$$ So, we can set $\alpha := ((p^{1/4}-1)(1+Ip^{1/4}))^{-1}$, which is a root of this equation. Then $(\star)$ turns into $$\left(\frac{i^2+j^2}p\right) = \frac{2}{I\sqrt{p}}\sum_{n=1}^{(p-1)/2} (\zeta^{i^2n^2}+\alpha)(\zeta^{j^2n^2}+\alpha),$$ implying that $$\left[\left(\frac{i^2+j^2}p\right)\right]_{i,j=1}^{(p-1)/2} = \frac{2}{I\sqrt{p}} A^2,$$ where matrix $A := \left[ \zeta^{i^2j^2}+\alpha \right]_{i,j=1}^{(p-1)/2}$. It follows that $$-S_p = T_p^2,\quad \text{where}\quad T_p := \det\left(\frac{1+I}{p^{1/4}}A\right).$$ The original question reduces to showing that $T_p$ is an integer.

This is not an answer, but a reduction to supposedly simpler problem. (edited with more details)

Using quadratic Gauss sums, we can express Legendre symbol $\left(\frac{i^2+j^2}p\right)$ as \begin{split} \left(\frac{i^2+j^2}p\right) &= \frac{1}{I\sqrt{p}}\sum_{n=0}^{p-1} \zeta^{(i^2+j^2)n^2} \\ &= \frac{1}{I\sqrt{p}}\left(1+2\sum_{n=1}^{(p-1)/2} \zeta^{i^2n^2}\zeta^{j^2n^2}\right),\qquad(\star) \end{split} where $\zeta := \exp(\frac{2\pi I}p)$ and $I$ is the imaginary unit. From $(\star)$ it follows that $$\left[\left(\frac{i^2+j^2}p\right)\right]_{i,j=1}^{(p-1)/2} = \frac{1}{I\sqrt{p}}(J + 2B^2),$$ where $J$ is the matrix of all ones and $B:=[\zeta^{i^2j^2}]_{i,j=1}^{(p-1)/2}$. This enables application of pavl0's approach in evaluating the determinant in order to show that it is a square.

Alternatively, $(\star)$ allows us to extract the square root from the determinant without computing its value. To do so, we need to embed "$1+$" into the quadratic sum in $(\star)$ by introducing a parameter $\alpha$: $$1+2\sum_{n=1}^{(p-1)/2} \zeta^{i^2n^2}\zeta^{j^2n^2} = 2\sum_{n=1}^{(p-1)/2} (\zeta^{i^2n^2}+\alpha)(\zeta^{j^2n^2}+\alpha).$$ Since $2\sum_{n=1}^{(p-1)/2} \zeta^{i^2n^2} = I\sqrt{p} - 1$ for any $i\not\equiv 0\pmod{p}$, we need $\alpha$ to satisfy the quadratic equation: $$1 = 2(I\sqrt{p} - 1)\alpha + (p-1)\alpha^2.$$ So, we can set $\alpha := ((p^{1/4}-1)(1+Ip^{1/4}))^{-1}$, which is a root of this equation. Then $(\star)$ turns into $$\left(\frac{i^2+j^2}p\right) = \frac{2}{I\sqrt{p}}\sum_{n=1}^{(p-1)/2} (\zeta^{i^2n^2}+\alpha)(\zeta^{j^2n^2}+\alpha),$$ implying that $$\left[\left(\frac{i^2+j^2}p\right)\right]_{i,j=1}^{(p-1)/2} = \frac{2}{I\sqrt{p}} A^2,$$ where matrix $A := \left[ \zeta^{i^2j^2}+\alpha \right]_{i,j=1}^{(p-1)/2}$. It follows that $$-S_p = T_p^2,\quad \text{where}\quad T_p := \det\left(\frac{1+I}{p^{1/4}}A\right).$$ The original question reduces to showing that $T_p$ is an integer.

This is not an answer, but a reduction to supposedly simpler problem.

Using quadratic Gauss sums, we can represent Legendre symbol $\left(\frac{i^2+j^2}p\right)$ as $$\left(\frac{i^2+j^2}p\right) = \frac{2}{I\sqrt{p}}\sum_{n=1}^{(p-1)/2} (\zeta^{i^2n^2}+a)(\zeta^{j^2n^2}+a),$$ where $\zeta := \exp(\frac{2\pi I}p)$, $a := ((p^{1/4}-1)(1+Ip^{1/4}))^{-1}$, and $I$ is the imaginary unit. That is, we have $$\left[\left(\frac{i^2+j^2}p\right)\right]_{i,j=1}^{(p-1)/2} = \frac{2}{I\sqrt{p}} A^2,$$ where matrix $A := \left[ \zeta^{i^2j^2}+a \right]_{i,j=1}^{(p-1)/2}$. It follows that $$-S_p = T_p^2,\quad \text{where}\quad T_p := \det\left(\frac{1+I}{p^{1/4}}A\right).$$ The original question reduces to showing that $T_p$ is an integer.

This is not an answer, but a reduction to supposedly simpler problem. (edited with more details)

Using quadratic Gauss sums, we can express Legendre symbol $\left(\frac{i^2+j^2}p\right)$ as \begin{split} \left(\frac{i^2+j^2}p\right) &= \frac{1}{I\sqrt{p}}\sum_{n=0}^{p-1} \zeta^{(i^2+j^2)n^2} \\ &= \frac{1}{I\sqrt{p}}\left(1+2\sum_{n=1}^{(p-1)/2} \zeta^{i^2n^2}\zeta^{j^2n^2}\right),\qquad(\star) \end{split} where $\zeta := \exp(\frac{2\pi I}p)$ and $I$ is the imaginary unit. From $(\star)$ it follows that $$\left[\left(\frac{i^2+j^2}p\right)\right]_{i,j=1}^{(p-1)/2} = \frac{1}{I\sqrt{p}}(J + 2B^2),$$ where $J$ is the matrix of all ones and $B:=[\zeta^{i^2j^2}]_{i,j=1}^{(p-1)/2}$. The enables application of pavl0's approach in evaluating the determinant in order to show that it is a square.

Alternatively, it is possible to "extract" the square root from the determinant without knowing its value. To do so, we need to embed "$1+$" into the quadratic sum in $(\star)$ by introducing a parameter $\alpha$: $$1+2\sum_{n=1}^{(p-1)/2} \zeta^{i^2n^2}\zeta^{j^2n^2} = 2\sum_{n=1}^{(p-1)/2} (\zeta^{i^2n^2}+\alpha)(\zeta^{j^2n^2}+\alpha).$$ Since $2\sum_{n=1}^{(p-1)/2} \zeta^{i^2n^2} = I\sqrt{p} - 1$ for any $i\not\equiv 0\pmod{p}$, we get that $\alpha$ satisfies the quadratic equation: $$1 = 2(I\sqrt{p} - 1)\alpha + (p-1)\alpha^2.$$ So, we can set $\alpha := ((p^{1/4}-1)(1+Ip^{1/4}))^{-1}$, which is a root of this equation. Then $(\star)$ turns into $$\left(\frac{i^2+j^2}p\right) = \frac{2}{I\sqrt{p}}\sum_{n=1}^{(p-1)/2} (\zeta^{i^2n^2}+\alpha)(\zeta^{j^2n^2}+\alpha),$$ implying that $$\left[\left(\frac{i^2+j^2}p\right)\right]_{i,j=1}^{(p-1)/2} = \frac{2}{I\sqrt{p}} A^2,$$ where matrix $A := \left[ \zeta^{i^2j^2}+\alpha \right]_{i,j=1}^{(p-1)/2}$. It follows that $$-S_p = T_p^2,\quad \text{where}\quad T_p := \det\left(\frac{1+I}{p^{1/4}}A\right).$$ The original question reduces to showing that $T_p$ is an integer.