If $p$ is a sum of two squares, let $p=x^2 + y^2$. Because $p$ must be odd, therefore one of $x$ and $y$ must be even and the other must be odd. Without loss of generality, let $x$ be even and $y$ be odd. Then it is easy to show that $x^2\equiv 0\pmod{4}$ and $y^2\equiv 1\pmod{4}$. Hence, $p=x^2 + y^2 \equiv 0 + 1\equiv 1\pmod{4}$.

On the other hand, if $p\equiv 1\pmod{4}$. We are going to show that $a^2 + 1=p$ is solvable for some integer $a$. If so, by # 4187, we know $p$ is a sum of two squares. For this, let's consider factors of $(p-1)!$: $$(p-1)!=1\cdot 2\cdot 3\cdots \left(\frac{p-1}{2}\right)\cdot\left(\frac{p+1}{2}\right)\cdots (p-2)\cdot (p-1)$$

Each factor in the first half can be paired with a corresponding one in the second half with negative congruent relationship: $$\begin{array}{rcll} 1&\equiv & -(p-1) &\pmod{p} \\ 2&\equiv & -(p-2) &\pmod{p} \\ & \cdots \\ \frac{p-1}{2}&\equiv & -\frac{p+1}{2} &\pmod{p} \end{array}$$

Therefore, we have $$(p-1)!\equiv (-1)^{\frac{p-1}{2}}\left(1\cdot 2\cdot 3\cdots \left(\frac{p-1}{2}\right)\right)^2\pmod{p}$$

Meanwhile, by Wilson's theorem, we have $(p-1)!\equiv -1\pmod{p}$. Thus $$-1 \equiv (-1)^{\frac{p - 1}{2}}\left(1\cdot 2\cdot 3\cdots \left(\frac{p-1}{2}\right)\right)^2\pmod{p}$$

or $$(-1)^{\frac{p + 1}{2}} \equiv \left(1\cdot 2\cdot 3\cdots \left(\frac{p-1}{2}\right)\right)^2\pmod{p}$$

As $p\equiv 1\pmod{4}$, let $p=4k+1$ where $k$ is an integer. Then $\frac{p+1}{2} = 2k+1$ is odd. Then the previous relation can be written as $$\left(1\cdot 2\cdot 3\cdots \left(\frac{p-1}{2}\right)\right)^2 + 1\equiv 0\pmod{p}$$

Let $a=1\cdot 2\cdot 3\cdots \left(\frac{p-1}{2}\right)$, we have $a^2 + 1\equiv 0\pmod{p}$ which is what we are looking for.