Reproducing OP's result
The notebook can be found on my github.
My version is "13.1.0 for Mac OS X ARM (64-bit) (June 16, 2022)"
and I can reproduce the OP's result by
Clear[f,I1]
f[x_,{r_,a_,L_}]:=(1-x^2)^L ((r+I a x)^(-5-L) (-((2+L) r)+2 I a x)+(-1)^(1+L) (r-I a x)^(-5-L) ((2+L) r+2 I a x))
I1[r_,a_,L_]:=Integrate[f[x,{r,a,L}],{x,-1,1},Assumptions->L>-1&&r>0&&a>0]
I1[r,a,2]
(*(64 (a^2-2 r^2))/(15 (a^2+r^2)^4)*)
I1[r,a,L]/.{L->2}//Simplify
(*-((8 (5 a^3+20 a^2 r+29 a r^2+16 r^3))/(15 r^5 (a+r)^4))*)
Actually if we modify the assumptions of I1
without $a>0$
I11[r_,a_,L_]:=Integrate[f[x,{r,a,L}],{x,-1,1},Assumptions->L>-1&&r>0]
Mathematica only returns results with $\Im a \neq 0$
I11[r,a,L]/.{L->2}//Simplify

Indeed there is an inconsistency of Mathematica.
A workaround
We need to absorb the factor I
into a
by
Clear[g,I2];
g[x_,{r_,a1_,L_}]:=Evaluate[f[x,{r,a,L}]/.a->-I a1]
I2[r_,a1_,L_]:=Integrate[g[x,{r,a1,L}],{x,-1,1},Assumptions->L>-1&&r>0]
Now the order of integration and limit is commutative:
I2[r,a1,2]/.a1->I a

I2[r,a1,L]
I2[r,a1,L]/.{L->2,a1->I a}//Simplify

Some explanations
As a multi-variable complex function, the hypergeometric ${}_p F_{q}(\{a\},\{b\},z)$ has branch cut at $z>1$ for generic values of $\{a\},\{b\}$, and has simple poles with respect to ${b}$. It's often hard to implement the analytic structures in symbolic computation.
The branch cut can be detected in the integral expression of the hypergeometric function as follows. For simplicity we only consider the ${}_2 F_{1}$ from dlmf
\begin{equation}
\mathbf{F}(a, b ; c ; z)=\frac{1}{\Gamma(b) \Gamma(c-b)} \int_0^1 \frac{t^{b-1}(1-t)^{c-b-1}}{(1-z t)^a} d t,
\end{equation}
The singularities at $t=0,1$ are okay for suitable choices of $b,c$. The only divergence happens at $t=1/z$. To see the branch cut we take $a=1$ and consider the distribution
\begin{equation}
\int^{1}_{0}dt f(t)\frac{1}{t-1/z}
\end{equation}
If $1/z\in (0,1)$ this is divergent except for test functions satisfying $f(1/z)=0$. But if $1/z$ has a small nonvanishing imaginary part $1/z -I \epsilon \in (0,1)$ it is convergent due to
\begin{equation}
\lim _{\epsilon \rightarrow 0^{+}} \frac{1}{x \pm i \epsilon}=\mp i \pi \delta(x)+\mathcal{P}\left(\frac{1}{x}\right)
\end{equation}
(This is called the regularization of distributions.)
The ambiguity of the sign of $\epsilon$ corresponds to the direction of $z$ approaching to the real interval, either from above or below, which means we have a branch cut of $f(z)$ at $z>1$.
Back to this question, I believe that starting from real axis and then analytic continue it onto the imaginary axis is more safe.
I1[L_, a_, r_] =
is changed toI1[L_, a_, r_] :=
, I get the same answers (theI2
answer). $\endgroup$