this is not super hard to prove. it is modeled on gauss(?) proof of fermats thm that primes of form 4k+1 are sums of two squares.
i.e. if p = 4k+1, then X^2+1 = 0 has a solution mod p. [e.g. mod 5, the solution is 1(2), and mod 13 the solution is 1(2)(3)(4)(5)(6) = 5, mod 13.]
this means that in (Z/p)[X], the polynomial X^2+1 is reducible, so that the ring (Z/p) is not a domain. but then in Z, the element p must be reducible, so p = (a+bi)(c+di), and taking absolute values of both sides, and squaring, we get p = a^2 + b^2.
Now Hurwitz defined the integral quaternions in a similar way, Z[i,j,k], (with a judicious denominator of 2), and it follows, similarly that no prime integer p is irreducible in there.
hence p = (a+bi+cj+dk)(e+fi+gj+hk), and again taking norms of both sides shows, with a little computation, that at least 4p is a sum of 4 squares.
but euler showed that if n is a sum mof 4 squares so is n/2, and applying this twice, p is a sum of 4 squares.
in both cases, the product of integers which are sums of 2 or 4 squares, is also. [this is because norms are multiplicative.] so the prime case does it all.
these proofs rely on the existence of a division algortihm in these two rings, even though the second one is non commutative.