Taking a look at the paper and the equations you mentioned in the comments, they're varying $I_{W_2}$ with respect to $g_{ab}$ rather than taking its covariant derivative. To see how you get Eq. 108 from Eq. 106, you either do the variation directly or integrate by parts in order for every term to be contracted with $\underline{\delta g^{ab}}$ while dismissing boundary terms. Especifically, you have
$$I_{W_2} = \int g^{ce} g^{df} R_{ef} R_{cd} \sqrt{-g} \text{d}^4 x.$$
Varying with respect to the metric gives
\begin{align}
\delta I_{W_2} &= \int \Big\{2R^c_{\ b} R_{ac} \sqrt{-g} \delta g^{ab} + 2 R^{cd} \delta R_{cd} \sqrt{-g} - \frac{1}{2} R_{cd} R^{cd} g_{ab} \sqrt{-g} \delta g^{ab} \Big\} \text{d}^4 x\\
&= \int \Big\{2R^c_{\ b} R_{ac} + 2 [\nabla_e (\delta \Gamma^e_{\ cd}) - \nabla_d (\delta \Gamma^e_{\ ec})] R^{cd} - \frac{1}{2} R_{cd} R^{cd} g_{ab} \delta g^{ab} \Big\} \sqrt{-g} \text{d}^4 x.
\end{align}
The variation of the Christoffel symbol is given by
\begin{align}
\delta \Gamma^e_{\ cd} &= \frac{1}{2} g^{ef} (\nabla_c \delta g_{df} + \nabla_d \delta g_{cf} - \nabla_f \delta g_{cd})\\
&= -\frac{1}{2} g^{ef} (g_{da} g_{fb} \nabla_c \delta g^{ab} + g_{ca} g_{fb} \nabla_d \delta g^{ab} - g_{ca} g_{db} \nabla_f \delta g^{ab}).
\end{align}
Focusing on $\nabla_e (\delta \Gamma^e_{\ cd}) R^{cd}$ gives
\begin{align}
\nabla_e (\delta \Gamma^e_{\ cd}) R^{cd} &= -\frac{1}{2} (R^c_{\ a} \nabla_b \nabla_c \delta g^{ab} + R_a^{\ \ d} \nabla_b \nabla_d \delta g^{ab} - R_{ab} \nabla^f \nabla_f \delta g^{ab})\\
&= -\frac{1}{2} (2 \nabla_c \nabla_b R^c_{\ a} - \square R_{ab}) \delta g^{ab} + \text{boundary terms}.
\end{align}
Similarly
$$\nabla_d (\delta \Gamma^e_{\ ec}) R^{cd} = -\frac{1}{2} (\nabla_b \nabla_c R^c_{\ \ a} + g_{ab} \nabla_c \nabla_d R^{cd} - \nabla_a \nabla_c R^c_{\ \ b}) \delta g^{ab} + \text{boundary terms}.$$
Plugging into $\delta I_{W_2}$ gives
$$\delta I_{W_2} = \int \Big\{2R^c_{\ b} R_{ac} + \square R_{ab} + \nabla_b \nabla_c R^c_{\ \ a} + g_{ab} \nabla_c \nabla_d R^{cd} - 2 \nabla_c \nabla_b R^c_{\ \ a} - \nabla_a \nabla_c R^c_{\ \ b} - \frac{1}{2} R_{cd} R^{cd} g_{ab} \Big\} \delta g^{ab} \sqrt{-g} \text{d}^4 x.$$
Using Bianchi's identity we obtain $\nabla_c \nabla_d R^{cd} = \frac{1}{2} \square R$, and $\nabla_b \nabla_c R^c_{\ \ a} - \nabla_c \nabla_b R^c_{\ \ a} = \frac{1}{2} (\nabla_b \nabla_a R - \nabla_a \nabla_b R) = 0$. Therefore
$$\delta I_{W_2} = \int \Big\{\frac{1}{2} \square R g_{ab} + \square R_{ab} - \nabla_c \nabla_b R^c_{\ \ a} - \nabla_a \nabla_c R^c_{\ \ b} + 2R^c_{\ b} R_{ac} - \frac{1}{2} R_{cd} R^{cd} g_{ab} \Big\} \delta g^{ab} \sqrt{-g} \text{d}^4 x,$$
which is Eq. 108. If there's some sign discrepancies is because I varied with respect to $g^{ab}$, not $g_{ab}$.
Edit: I realize that, strictly speaking, I didn't answer your question. However, here's what I would do. As you correctly point out, $\sqrt{-g}$ is a tensor density, which transforms as
$$\nabla_a \sqrt{-g} = \partial_a \sqrt{-g} - \Gamma^b_{\ ab} \sqrt{-g}.$$
Now, given a basis $\{ \partial_\mu \}$ and using the identity
$$\ln (-g) = \text{Tr} (\ln ([g_{\mu \nu}]))$$
we obtain
$$\partial_\alpha \ln (-g) = g^{\mu \nu} \partial_\alpha g_{\mu \nu}.$$
It's direct to see $\Gamma^\beta_{\ \alpha \beta} = \frac{1}{2} g^{\beta \gamma} \partial_\alpha g_{\beta \gamma}$. This implies $\Gamma^\beta_{\ \alpha \beta} = \frac{1}{2} \partial_\alpha \ln (-g) = \partial_\alpha \ln (\sqrt{-g}).$ Plugging the former in $\nabla_a \sqrt{-g}$ results in
\begin{align}
\nabla_\alpha \sqrt{-g} &= \partial_\alpha \sqrt{-g} - \partial_\alpha \ln (\sqrt{-g}) \cdot \sqrt{-g}\\
&= \partial_\alpha \sqrt{-g} - \frac{\partial_\alpha \sqrt{-g}}{\sqrt{-g}} \sqrt{-g}\\
&= 0.
\end{align}
I'm sure the paper's author meant varying with respect to the metric instead of taking the covariant derivative in Eq. 106. In fact, that is what they're doing. Otherwise, it doesn't make sense to take the covariant derivative of this equation either within the context of the paper or the author's reply to your email.
