Explicit formula for Fermat's 4k+1 theorem
Here is a high level proof. I assume it can be done in a more elementary way. Chapter 3 of Silverman's Arithmetic of Elliptic Curves is a good reference for the ideas I am using.
Let $E$ be the elliptic curve $y^2 = x^3+x$. By a theorem of Weyl, the number of points on $E$ over $\mathbb{F}_p$ is $p- \alpha- \overline{\alpha}$ where $\alpha$ is an algebraic integer satisfying $\alpha \overline{\alpha} =p$, and the bar is complex conjugation. (If you count the point at $\infty$, then the formula should be $p - \alpha - \overline{\alpha} +1$.)
Let $p \equiv 1 \mod 4$. We will establish two key claims: Claim 1: $\alpha$ is of the form $a+bi$, for integers $a$ and $b$, and Claim 2: $-2a \equiv \binom{(p-1)/2}{(p-1)/4} \mod p$. So $a^2+b^2 = p$ and $a \equiv -\frac{1}{2} \binom{(p-1)/2}{(p-1)/4}$, as desired.
Proof sketch of Claim 1: Let $R$ be the endomorphism ring of $E$ over $\mathbb{F}_p$. Let $j$ be a square root of $-1$ in $\mathbb{F}_p$. Two of the elements of $R$ are $F: (x,y) \mapsto (x^p, y^p)$ and $J: (x,y) \mapsto (-x,jy)$.
Note that $F$ and $J$ commute; this uses that $j^p = j$, which is true because $p \equiv 1 \mod 4$. So $F$ and $J$ generate a commutative subring of $R$. If you look at the list of possible endomorphism rings of elliptic curves, you'll see that such a subring must be of rank $\leq 2$, and $J$ already generates a subring of rank $2$. (See section 3.3 in Silverman.) So $F$ is integral over the subring generated by $J$. That ring is $\mathbb{Z}[J]/\langle J^2=-1 \rangle$, which is integrally closed. So $F$ is in that ring, meaning $F = a+bJ$ for some integers $a$ and $b$.
If you understand the connection between Frobenius actions and points of $E$ over $\mathbb{F}_p$, this shows that $\alpha = a+bi$.
Proof sketch of Claim 2: The number of points on $E$ over $\mathbb{F}_p$ is congruent modulo $p$ to the coefficient of $x^{p-1}$ in $(x^3+x)^{(p-1)/2}$ (section 3.4 in Silverman). This coefficient is $\binom{(p-1)/2}{(p-1)/4}$. So $$- \alpha - \overline{\alpha} \equiv \binom{(p-1)/2}{(p-1)/4} \mod p$$ or $$-2a \equiv \binom{(p-1)/2}{(p-1)/4} \mod p$$ as desired.
Remark: This is very related to the formula Matt E mentions. For $u \in \mathbb{F}_p$, the number of square roots of $u$ in $\mathbb{F}_p$ is $1+\left( \frac{u}{p} \right)$. So the number of points on $E$ is $$p+\sum_{x \in \mathbb{F}_p} \left( \frac{x^3+x}{p} \right).$$
This is essentially Matt's sum; if you want, you could use the elliptic curve $y^2 = x^3-x$ in order to make things exactly match, although that would introduce some signs in other places. So your remark gives another (morally, the same) proof of Claim 2.
There is a proof on page 192 of Franz Lemmermeyer's book, Reciprocity Laws from Euler to Eisenstein. I found it by typing "sum of two squares" and "binomial coefficient" into Google.
There is also a proof in Allan Adler's paper, Eisenstein and the Jacobian varieties of Fermat curves, which paper is freely available on the web.
Denote by $\phi(D)$ the Jacobsthal sum $$ \sum_t \left(\frac tp\right)\left(\frac{t^2+D}p\right). $$ Note that $\phi(a^2D)=(a/p)\phi(D)$, so $|\phi(D)|$ depends only on $\bigl(\frac Dp\bigr)$; let $2x$ and $2y$ be the absolute values of the Jacobsthal sum for quadratic residue and quadratic non-residue respectively.
Theorem (Jacobsthal, 1907). If $p=4k+1$, $x^2+y^2=p$.
(As explained in the OP, this statement implies the formula of Gauss.)
Proof. Obviously, $$ \sum_{D\in\mathbb F_p}\phi(D)^2=2(p-1)(x^2+y^2), $$ so we need to compute the sum \begin{align*} \sum_{D}\phi(D)^2 & = \sum_{D,s,t}\left(\frac{st}p\right)\left(\frac{s^2+D}p\right)\left(\frac{t^2+D}p\right) \\ & = \sum_{s,t}\left(\frac{st}p\right)\sum_D\left(\frac{s^2+D}p\right)\left(\frac{t^2+D}p\right). \end{align*}
Lemma. $$ \sum_D\left(\frac Dp\right)\left(\frac{D+c}p\right)=\begin{cases} \phantom p-1,&c\neq 0;\\ p-1,&c=0. \end{cases} $$
To prove the lemma for $c\neq0$ observe that in this case \begin{align*} \sum_D\left(\frac{D^2+cD}p\right) &= \sum_{D\neq0}\left(\frac{1+cD^{-1}}p\right)= \sum_{D'\neq1}\left(\frac{D'}p\right) \\ & =-\left(\frac1p\right)=-1. \end{align*}
Now let us apply the lemma to our sum: $$ \sum_D\left(\frac{s^2+D}p\right)\left(\frac{t^2+D}p\right)= \sum_{D'}\left(\frac{D'}p\right)\left(\frac{D'+t^2-s^2}p\right)= \begin{cases} \phantom p-1,&s\neq\pm t;\\ p-1,&s=\pm t; \end{cases} $$ so $$ \sum_{D}\phi(D)^2=(p-1)\left\{\sum_{s=t}\left(\frac{t^2}p\right)+\sum_{s=-t}\left(\frac{-t^2}p\right)\right\}+(-1)\sum_{s\neq\pm t}\left(\frac{st}p\right). $$ The sums in the brackets $\{ \}$ clearly are equal to $2(p-1)$; the last sum is $-2(p-1)$ because $\sum (\frac{st}p) = 0$, so $$ \sum_{D}\phi(D)^2 = 2(p-1)^2 + 2(p-1) = 2(p-1)p. $$ Hence indeed $x^2+y^2=p$.
(To my surprise, although the proof is quite short and elementary, I haven't been able to find a concise modern reference. Anyway, Jacobsthal's paper «Über die Darstellung der Primzahlen der Form $4n+1$ als Summe zweier Quadrate» is quite accessible.)