We can use the excellent idea of the answer here.
The idea is to show that the second derivative is negative, which implies the function is concave.
Define $\widetilde{A}=A+\epsilon_1 C$ and $\widetilde{B}=B+\epsilon_2 D$, with $C,D$ symmetric.
So, we should show:
$$
\frac{\partial^2}{\partial \epsilon_i^2}
\text{tr}
\left( \left[\widetilde{A}^{-1} +\widetilde{B}^{-1}\right]^{-1}\right)\bigg|_{\epsilon_i=0}
\leq 0
$$
Notice that we can do it for $i=1$, then the result holds by symmetry.
Next we will use the Woodbury matrix formula and the Taylor expansion of the inverse function for matrices (see here as well).
Note the assumption that $||A^{-1}B||<1$.
This lets us write:
\begin{align*}
\left[ \widetilde{A}^{-1} +\widetilde{B}^{-1} \right]^{-1}
&=
\widetilde{A} - \widetilde{A}(\widetilde{A} + \widetilde{B})^{-1}\widetilde{A}\\
&=
\widetilde{A} -
\widetilde{A}
(\widetilde{A}^{-1} - \widetilde{A}^{-1}\widetilde{B}\widetilde{A}^{-1}+
\widetilde{A}^{-1}\widetilde{B}\widetilde{A}^{-1}\widetilde{B}\widetilde{A}^{-1}-\ldots\,)
\widetilde{A}\\
&=
\widetilde{B} - \widetilde{B}\widetilde{A}^{-1}\widetilde{B}+
\widetilde{B}\widetilde{A}^{-1}\widetilde{B}\widetilde{A}^{-1}\widetilde{B}+\ldots\\
&=
\widetilde{B}
-
\widetilde{B}
\left[
{A}^{-1} - \epsilon_1{A}^{-1}C{A}^{-1}+
\epsilon^2_1{A}^{-1}C{A}^{-1}C{A}^{-1}+
\ldots
\right]
\widetilde{B}
+ \ldots
\end{align*}
Thus we can get:
\begin{align*}
\frac{\partial^2}{\partial \epsilon_i^2}
\text{tr}
\left( \left[\widetilde{A}^{-1} +\widetilde{B}^{-1}\right]^{-1}\right)\bigg|_{\epsilon_i=0}
&=
\text{tr}(-2\widetilde{B}
\underbrace{
A^{-1}CA^{-1}CA^{-1}}_{PA^{-1}P^T}
\widetilde{B})\\
&=-2 \text{tr}(-2\widetilde{B}
\underbrace{PA^{-1}P^T}_{\text{Pos. Semi. Def}}
\widetilde{B})\\
&\leq 0
\end{align*}
using the fact that $B$ and $A^{-1}$ are positive definite.
Since the second derivative is non-positive, the trace function
$\text{tr}([A^{-1}+B^{-1}]^{-1})$
is concave with respect to a perturbation of its arguments.