How can one prove that numbers of the form $4x^2+1$ can only be divided by primes of the form $4y+1$ (e.g. there is no $x$ for which $7$ divides $4x^2+1$)?
On a quick lookup, the statement is given as the first comment in https://oeis.org/A053755 and of course, manual testing concurs with this
In general; note that your question is about the form $4x^2 + y^2$ with $\Delta = -16.$ If $q \equiv 3 \pmod 4$ is a prime that divides the form, then it divides both $x,y.$ However, in your case $y=1.$
Given a binary quadratic form $$ f(x,y) = a x^2 + b xy+ c y^2 $$ with $a,b,c$ integers. Given $$ \Delta = b^2 - 4 a c, $$ where we require that $\Delta,$ if non-negative, is not a square (so also $\Delta \neq 0,1$).
Proposition: given an odd prime $q$ such that $q$ does not divide $\Delta$ and, in fact, $$ (\Delta| q) = -1, $$ whenever $$ f(x,y) \equiv 0 \pmod q, $$ then BOTH $$ x \equiv 0 \pmod q, \; \; \; \; \; y \equiv 0 \pmod q. $$
Proof: the integers taken $\pmod q$ form a field; every nonzero element has a multiplicative inverse. If either $a,c$ were divisible by $q,$ we would have $\Delta $ equivalent to $b^2$ mod $q,$ which would cause $(\Delta|q)$ to be $1$ rather than $-1.$ As a result, neither $a$ nor $c$ is divisible by $q.$
Next, $q$ is odd, so that $2$ and $4$ have multiplicative inverses mod $q,$ they are non zero in the field. Put togethaer, $$4a \neq 0 \pmod q$$
Now, complete the square: $$ a x^2 + b xy+ c y^2 \equiv 0 \pmod q $$ if and only if $$ 4a(a x^2 + b xy+ c y^2) \equiv 0 \pmod q, $$ $$ 4a^2 x^2 + 4abxy + 4ac y^2 \equiv 0 \pmod q, $$ $$ 4a^2 x^2 + 4abxy + (b^2 y^2 - b^2 y^2) + 4ac y^2 \equiv 0 \pmod q, $$ $$ 4a^2 x^2 + 4abxy + b^2 y^2 - (b^2 y^2 - 4ac y^2) \equiv 0 \pmod q, $$ $$ (4a^2 x^2 + 4abxy + b^2 y^2) - (b^2 - 4ac) y^2 \equiv 0 \pmod q, $$ $$ (2ax + by)^2 - (b^2 - 4ac) y^2 \equiv 0 \pmod q, $$ $$ (2ax + by)^2 - \Delta y^2 \equiv 0 \pmod q, $$ $$ (2ax + by)^2 \equiv \Delta y^2 \pmod q. $$
Now, ASSUME that $y \neq 0 \pmod q.$ Then $y$ has a multiplicative inverse which we are allowed to call $1/y,$ and we have $$ \frac{(2ax + by)^2}{y^2} \equiv \Delta \pmod q, $$ $$ \left( \frac{2ax + by}{y} \right)^2 \equiv \Delta \pmod q. $$ However, the HYPOTHESIS that $ (\Delta| q) = -1 $ says that this is impossible, thus contradicting the assumption that $y \neq 0 \pmod q.$
So, in fact, $y \equiv 0 \pmod q.$ The original equation now reads $$ a x^2 \equiv 0 \pmod q $$ with the knowledge that $a \neq 0 \pmod q,$ so we also get $$ x \equiv 0 \pmod q, \; \; \; \; \; y \equiv 0 \pmod q. $$ $$ \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc \bigcirc $$