As far as I understand mathematical intuition, writing down such a formula without "walking" through the usual technical road to it must be the result of some idea, especially in mathematical areas that are considered frontiers (and quaternions were indeed a a frontier of mathematics at the times of Gauss).
Example: My unsuccesful interpretation
I mention this unsuccesful attempt because maybe someone will be able to continue this line of thought and rewrite it in the language of modern abtract algebra.
The proof above views quaternions as ordered pairs of complex numbers (for example: $q_1 = x+jy$ where $x,y$ are complex numbers); Gauss has already encountered the phenomenon of quaternionic behaviour in ordered pairs of complex numbers in another place (Gauss's werke, vol 3, p.384). Therefore I have an intuitive feeling he casted his statements in this particular form because he viewed the quaternions as an hypercomplex number system over the complex numbers (as an extension of $\mathbb{C}$), in the same way complex numbers can be viewed as an hyperreal number system over the reals.
Therefore I thought it is logical to try and search for a similar pattern in the complex numbers. Let $$c_1=a+bi,c_2=\alpha+\beta i$$ $$c_3 = c_1\cdot c_2 =A+Bi = (a\alpha - b\beta)+(b\alpha+a\beta)i$$
and then one gets that $$\frac{B}{A}-\frac{b}{a} = \frac{(b\alpha+a\beta)a-(a\alpha - b\beta)b}{Aa} = \frac{(a^2+b^2)\beta}{Aa} \equiv 0 \pmod {a^2+b^2}$$
where the last relation holds for prime values of $a^2+b^2$. The fact that this congruence is also correct for complex numbers is obvious from the proof for general quaternions, since the complex integers are a subset of quaternion integers. But maybe by looking at Gauss's congruence in the simpler case of $\mathbb{C}$ one can get idea about the kind of notions Gauss attempted to generalize from $\mathbb{C}$ to the quaternions.