What effect does rotation in the spatial domain has on phase in Fourier transforms?
This question has a simpler answer for the 2-D continuous-space Fourier transform but itsdDiscrete Fourier transform based verification requires some elaboration and careful implementation as @MarcusMüller has already mentioned.

On the continuous-space it can be shown that rotating an image $f(x,y)$ by $\theta$ radians CCW (counter-clock wise) on the $xy$ plane, will also rotate the corresponding CSFT (continuous-space Fourier transform) $F( \Omega_1 , \Omega_2 )$ by $\theta$ radians CCW :

CCW rotation by $\theta$ on the $xy$ plane can be described by the following transform:

$$
x' = x \cos(\theta) - y \sin(\theta) \\
y' = x \sin(\theta) + y \cos(\theta) \tag{1}
$$ 

where $x,y$ are the original coordinates and $x',y'$ are the rotated (new) coordinates. A corresponding inverse rotation by $\theta$, (or CCW by $-\theta$ ) has the following transform :

$$
x = + x' \cos(\theta) + y' \sin(\theta)  \tag{2} \\
y = - x' \sin(\theta) + y' \cos(\theta) 
$$ 

Given the original image $f(x,y)$, we rotate it CCW by $\theta$ to obtain the rotated image $g_r(x,y)$, then the following will hold :

$$ g_r(x,y) = f( x \cos(\theta) + y \sin(\theta) , - x \sin(\theta) + y \cos(\theta) ) \tag{3} $$ 

Note that in Eq. (3), the arguments of the original function $f(\cdot,\cdot)$ are described by the inverse rotation relation to coordinates $x,y$ of $g_r(\cdot,\cdot)$.

The CSFT of the rotated image $g_r(x,y)$ is:

$$
\boxed{ G_r(\Omega_1,\Omega_2) = \iint_{-\infty}^{\infty} g_r(x,y) e^{-j( \Omega_1 x + \Omega_2 y) } dx dy} \\
$$

$$
\begin{align} 
&amp;= \iint f( x \cos(\theta) + y \sin(\theta) , - x \sin(\theta) + y \cos(\theta) ) e^{-j( \Omega_1 x + \Omega_2 y) }   \tag{4} \\
&amp;= \iint f( x' , y' ) e^{-j[ \Omega_1 (x' \cos(\theta) - y'\sin(\theta)) + \Omega_2 (x' \sin(\theta) + y' \cos(\theta)) ]} dx' dy'  \tag{5}\\
&amp;= \iint f( x' , y' ) e^{-j[ ( \Omega_1 \cos(\theta) + \Omega_2  \sin(\theta) ) x' + (-\Omega_1 \sin(\theta) + \Omega_2 \cos(\theta)) y' ]} dx' dy'  \tag{6}\\
&amp;= F(  \Omega_1 \cos(\theta) + \Omega_2  \sin(\theta) ,   -\Omega_1 \sin(\theta) + \Omega_2 \cos(\theta) ) \tag{7}
\end{align}
$$

In moving from step (4) to (5), make the substitution defined in Eqs (1) and (2), and note that the Jacobian of the transformation is unity. See page 778 of Calculus Adams 6E for an explanation and derivation of this result.

In moving from step (6) to (7) we simply recognise the expression in Eq.(6) as the CSFT of the original signal $f(x,y)$ with arguments given according to relation in Eq.(3) given as below:

$$G_r(\Omega_1,\Omega_2) = F(  \Omega_1 \cos(\theta) + \Omega_2  \sin(\theta) ,   -\Omega_1 \sin(\theta) + \Omega_2 \cos(\theta) ) \tag{8} $$

Eq.(8) is analogous to Eq.(3) and states that the CSFT of the rotated image is also roated by the same amount;  $\theta$ CCW. This, therefore, means that phase is also rotated by the same amount.

Verification of this result within the discrete-space sequences $f[n,m]$ and their corresponding DFT sequences requires specific care given to potential aliasing caused by inadequate sampling of the rotated image $g_r(x,y)$, and  proper circular extension of the rotated sequence $g_r[n,m]$. 

Without going into mathematical details, I'll put the results of the following implementation of the rotation using Matlab/ Octave (rotation is wrt upper left corner)



We can still see artefacts of non-ideal implementation...

Note that with typical images it's harder to realize the ideal result.