Wilson lines indeed act as sources for charged particles. For the purposes of illustration, let me go to the abelian case momentarily, so that $\mathrm{tr}\,\mathrm{P}\!\exp(\cdot)$ is simply $\exp(\cdot)$, and I'll return to Yang-Mills in the end.
$\newcommand{\d}{\mathrm{d}}$
Abelian case
So consider an abelian gauge theory with a gauge field, $a$, and now you want to add a source. The obvious generalisation of "inserting some $\exp(\mathrm{i}\;J\cdot\phi)$ in the path integral" is inserting some
$$\exp\!\left(\mathrm{i}\;\int_\mathcal{M}\d^d x\ J^\mu(x)\ a_\mu(x)\right)\tag{1}$$
in the path integral. And now, let's come to differential forms. The gauge field, $a = a_\mu(x)\;\d x^\mu$, is a one-form and so is $J = J_\mu(x)\; \d x^\mu$. The integral in the expression (1) can be rewritten as
$$ \int_\mathcal{M}\d^d x\ J^\mu(x)\ a_\mu(x) = \int_\mathcal{M} a\wedge\star J =: \int_\mathcal{M} a\wedge K,$$
where $\star$ is the Hodge-star operator and I defined the $(d-1)$-form
$$K :=\star J = K_{\mu_1 \mu_2\cdots \mu_{d-1}} \d{x}^{\mu_1} \wedge \d{x}^{\mu_2} \wedge \cdots \wedge \d{x}^{\mu_{d-1}}.$$
Now we come to the fun part, Poincaré duality! Since your theory has gauge invariance, $a\sim a+\d \lambda $ an integration by parts shows that $\d\star J = 0$, or equivalently $\d K = 0$. That means that $K$ is a representative of a cohomology class of the $(d-1)$-th cohomology group of $\mathcal{M}$. Invoking Poincaré duality, i.e. $\mathrm{H}^n(\mathcal{M}) \cong \mathrm{H}_{d - n}(\mathcal{M})$ you can rewrite the above integral as
$$ \int_\mathcal{M} a\wedge K = \int_\widehat{K} a, \tag{2}$$
where $\widehat{K}$ is the Poincaré dual 1-cycle to the $(d-1)$-class $[K]$ (see appendix for a proof of (2)). That's it! If I rename $\widehat{K}$ as $C$, this argument just shows that
$$\exp\!\left(\mathrm{i} \int_C a\right) =: W[C],$$
is a source for charged particles.
In the opposite direction, starting with a Wilson loop $W[C]$ you can reverse engineer $J$ from (1), as
$$J_\mu(x) = (-1)^{(d-1)+\epsilon}\big[\star\!\widehat{C}\,\big]_\mu(x),$$
where $\epsilon=0$ if you are in Euclidean signature and $\epsilon=1$ in Lorentzian, and $\widehat{C}$ is (a representative of) the Poincaré dual $(d-1)$-class, of the 1-cycle $C$. This is roughly speaking a $(d-1)$-form delta function, restricting you to lie in $C:\quad$ $ \widehat{C} = \delta_{[d-1]}\big(x\in C\big), $ so
$$ J_\mu(x) = (-1)^{(d-1)+\epsilon} \left[\star\delta_{[d-1]}\big(x\in C\big)\right]_\mu.$$
Back to Yang-Mills
All of the above differential-geometric arguments hold the same in a non-abelian gauge theory, by replacing exponentials with traces of path-ordered exponentials, therefore the slogan is inserting a Wilson loop
$$W[C]:= \mathrm{tr}\,\mathrm{P}\!\exp\!\left(\mathrm{i}\int_C a\right)$$
in the path-integral, is equivalent to adding a source given by a current:
$ J_\mu(x) = (-1)^{(d-1)+\epsilon} \left[\star\delta_{[d-1]}\big(x\in C\big)\right]_\mu.$ If you wish to find the ground state energy in the presence of a Wilson loop, everything is now ready to be translated into Hamiltonian language, exactly like you did for the scalar case.
The path-ordering symbol somewhat obscures the source interpretation. One way to get around this, is described also in Tong's gauge theory lecture notes, around subsection 2.1.3. The trick is to trade the path-ordering symbol for auxiliary fields. In particular, you can write:
\begin{align} W_\mathrm{R}[C] &= \mathrm{tr}_\mathrm{R}\,\mathrm{P}\!\exp\!\left(\mathrm{i}\int_C a\right) \\
&= \int \mathrm{D}\lambda \mathrm{D}w^\dagger \mathrm{D}w \ \mathcal{O}_\mathrm{R}\left(w,w^\dagger\right)\ \exp\!\left(\mathrm{i}\int_C \left(w^\dagger \mathrm{i}\,\d w + \lambda (w^\dagger w-\kappa_\mathrm{R}) + w^\dagger w\; a\right) \right), \tag{3}
\end{align}
where $w$ is a complex vector of dimension $\dim\mathrm{R}$ ($\mathrm{R}$ is the representation on which the Wilson loop is computed), $\kappa_\mathrm{R}$ is number that depends on the representation, that controls the length of $w$, and $\mathcal{O}_\mathrm{R}$ is an insertion whose exact form depends on the representation. From here on, you can play the Poincaré game, make $C$ into a $J$ and you are left with a source term [the last term in (3)] and a bunch of auxiliary fields that you need to integrate out.
Appendix: Justifying Poincaré duality
$\newcommand{\cM}{\mathcal{M}}$In this appendix I will justify equation (2), namely, I will show that Poincaré duality can be used with $a$ being any $1$-form. Let's define the pairing
\begin{align}
\Omega^1(\cM)\times\Omega^{d-1}(\cM)&\to \mathbb{R} \\
\left(\omega,\eta\right)&\mapsto \int_\cM \omega\wedge\eta.
\end{align}
Fixing an element $a\in\Omega^1(\cM)$ and restricting the second factor to a cohomology class, induces a map
\begin{align}
\phi_a:\mathrm{H}^{d-1}(\cM) &\to \mathbb{R} \\
\eta &\mapsto \phi_a(\eta) = \int_\cM a\wedge \eta.
\end{align}
I can also construct a map
\begin{align}\psi_a :\mathrm{H}_{1}(\cM) &\to \mathbb{R} \\
C&\mapsto \psi_a(C):=\int_C a,
\end{align}
which is simply integrating $a$ on a 1-cycle. Poincaré duality asserts that $$\mathrm{H}_{1}(\cM) \cong \mathrm{H}^{d-1}(\cM).$$
I can rewrite that by saying that there is an isomorphism
$$\xi :\mathrm{H}_{1}(\cM) \overset{\sim}{\longrightarrow} \mathrm{H}^{d-1}(\cM).$$
Finally I can construct the following commutative square of maps:
$$\require{AMScd}
\begin{CD}
\mathrm{H}_{1}(\cM) @>{\xi}>{\sim}> \mathrm{H}^{d-1}(\cM) \\
@V{\psi_a}VV @VV{\phi_a}V \\
\mathbb{R}@>{\sim}>{\mathrm{id}}>\mathbb{R},
\end{CD}$$
which shows that $\phi_a\circ\xi = \mathrm{id}\circ\psi_a=\psi_a$, or in other words, for any $a\in\Omega^1(\cM)$
$$\int_C a=:\psi_a(C) = \phi_a(\xi(C)) := \int_\cM a\wedge \xi(C).$$
This is precisely equation (2) with $\xi(\widehat{K})=K$.
