The missing conceptual point to appreciate is that when you are calculating a correlation function it can be useful to interpret the same expression in different ways.
We wish to calculate the correlation function with a path integral by
$\left\langle U_g(\Sigma) W(q,\gamma) \right\rangle = \int \mathcal{D}(\text{fields}) U_g(\Sigma) W(q,\gamma) e^{iS}$
where $S$ is the action of the theory. We can interpret this as an expectation value of $U_g(\Sigma) W(q,\gamma)$ in the theory defined by $S$, so in the path integral, the contributions of different paths are weighted by $e^{iS}$. Maybe it's useful to keep track of that $\langle U_g(\Sigma) W(q,\gamma) \rangle_S$.
But since the Wilson line operator is itself $W(q,\gamma) = e^{i q \int_\gamma A}$, it can be convenient to simply reinterpret this term as part of the action,
$\int \mathcal{D}(\text{fields}) U_g(\Sigma) e^{i q \int_\gamma A} e^{iS} = \int \mathcal{D}(\text{fields}) U_g(\Sigma) e^{iS+i q \int_\gamma A}$.
Of course this equality is trivial. The only further trick used is to write this as an integral over a local Lagrangian density by making use of a delta function localized on $\gamma$, $\int_\gamma A = \int_{M^4} \delta^{(3)}(\gamma) A$.
So we've shown then that we can alternatively interpret this correlator in the theory $S$ as a correlator in the theory defined by $S'$,
$\langle U_g(\Sigma) W(q,\gamma) \rangle_S = \langle U_g(\Sigma) \rangle_{S'}$,
where if $S = \int_{M^4} \mathcal{L}$, then $S' = \int_{M^4} (\mathcal{L} + \delta^{(3)}(\gamma) A)$.
This reinterpretation is then useful because we can evaluate the path integral easily now by looking at the modified equations of motion in the theory $S'$, where taking the variation $\delta/\delta A$ now yields an extra term in the equation of motion.
Inserting that modified equation of motion into the operator $U_g(\Sigma)$ then gives us the result for that correlator. We interpret $\oint_{\Sigma_2} \star F$ via Stokes' theorem as $\int_{\Sigma_3} d \star F$ for a 3-surface such that $\partial \Sigma_3 = \Sigma_2$, and $\int_{\Sigma_3} \delta^{(3)}(\gamma)$ counts the linking of $\gamma$ and $\Sigma_2$. Going through the same steps with $U_g(\Sigma'_2)$ where $\Sigma'_2$ by definition does not link $\gamma$ then gives precisely their Eqn 2.43.
If it would be useful to see just such a calculation done with more explicit indices and a more familiar physical interpretation see Appendix B.5 of Koren (2022). The interpretation here of the calculation is the Aharonov-Bohm phase of a Wilson line around a cosmic string in $\mathbb{Z}_N$ gauge theory, $\langle e^{i\int_\gamma A}e^{i \int_{\Sigma_2} B} \rangle$. Here $B$ is the two-form gauge field which appears in the BF description of $\mathbb{Z}_N$ gauge theory. Going through the calculation there you can then clearly see that moving $ e^{i \int_{\Sigma_2} B} $ into the action introduces a source of magnetic flux, and further evaluating $e^{i\int_\gamma A}$ in that background gives the Aharonov-Bohm phase.
Actually the two computations are more closely related than you might think, since as discussed in Sec 3.2.3 of Brennan & Hong, in the discrete case the Wilson loop of the two-form gauge symmmetry is itself really the symmetry defect operator for the one-form global symmetry. Thus because of this interesting duality of discrete gauge theories, this computation has exactly the sort of interpretation as the one in the question. Finally, let me also mention that Cheung, Derda, Kim, Nevoa, Rothstein, Shah (2024) also has some nice pedagogical discussions in this direction.
