OP Here, I've figured out the issue. It really is a functional derivative. Below I include an excerpt from my own notes, in case another lost soul stumbles upon this page in the future. (Please excuse the changes in notation from my original post, I hope it won't be difficult to switch over, and please excuse the fact that some of this might be in excruciating detail - also certainly not the most efficient way of deriving the result)
$\newcommand{\D}{\partial}$
$\newcommand{\dd}{d}$
We've already encountered Wilson lines through a Feynman diagram analysis. It will turn out that they are a common object that arises when computing the propagator in the presence of a background field. In the most general form, Wilson lines are:
\begin{align}
V(z,y) \equiv \hat{P} \exp \left(ig \int_{y}^{z} A_\mu(s)\dd s^\mu\right).
\end{align}
$\hat{P}$ is the path ordering operator. The integral here is a line integral in spacetime from some starting point $y^\mu$ to some ending point $z^\mu$, along some path that we can parametrize via $s^\mu(t)$.
\begin{align}
V(z,y) = \hat{P} \exp \left( ig \int_0^1 A_\mu(s(t)) \frac{ \dd s^\mu}{\dd t} \dd t\right).
\end{align}
We've defined $s^\mu(0) = y^\mu$ and $s^\mu(1) = z^\mu$. In order to see how this might appear as a solution to QCD equations of motion, we'll need to know how to take derivatives of Wilson lines. We assume here that the background field $A_\mu$ is a fixed function, and therefore the Wilson line is a functional of the path $s^\mu(t)$. We want to take a derivative with respect to spacetime coordinates $\frac{ \D }{\D x^\mu}$. Since $V(z,y)$ is a functional, we can think of it as a function of an infinite number of variables $s_i^\mu$. The label $i$ is purely symbolic, since it's really a function of $s(t)$ for every value of $t$. Using this symbolic notation, one can then use the chain rule:
\begin{align}
\frac{ \D V(z,y)}{\D x^\mu} = \frac{ \D V(z,y)}{\D s_i^\nu}\frac{ \D s_i^\nu}{\D x^\mu}.
\end{align}
Moving back to the functional notation, we can write $s_i^\nu = s^\nu(t)$ and then the implied sum over $i$ becomes an integral over $t$.
\begin{align} \label{Dform}
\frac{ \D V(z,y)}{\D x^\mu} = \int_0^1 \dd t \frac{ \delta V(z,y)}{\delta s^\nu(t)} \frac{ \D s^\nu}{\D x^\mu}.
\end{align}
Now it is clear that we cannot go further without figuring out the functional derivative of the Wilson line with respect to the path. Let's recall the method of finding the functional derivative. First, let's write $V(z,y)$ as a functional of $s^\nu(t)$.
\begin{align}
V[s(t)] = \hat{P} \exp\left( ig\int_0^1 \dd t A_\mu(s(t))\frac{ \dd s^\mu}{\dd t}\right).
\end{align}
Then let's make the transformation $s^\mu(t) \rightarrow s^\nu(t) + \delta s^\mu(t)$. We'll define $\delta V = V(s+\delta s) - V(s)$. Then, all we have to do is find $\delta V$ by simplifiying and ommitting any factors of $\mathcal{O}(\delta^2)$ or higher, and write the result as:
\begin{align} \label{funcderiv}
\delta V = \int_0^1 \dd t \frac{ \delta V}{\delta s^\mu(t)} \delta s^\mu(t),
\end{align}
and this will suffice to identify $ \frac{ \delta V}{\delta s^\mu(t)}$. So, let's begin as follows:
\begin{align}
\delta V = \hat{P}\exp\left(ig \int_0^1 \dd t A_\mu(s+\delta s)\left( \frac{ \dd s^\mu}{\dd t} + \frac{ \dd \delta s^\mu}{\dd t}\right)\right) - \hat{P} \exp\left(ig\int_0^1\dd t A_\mu(s) \frac{ \dd s^\mu}{\dd t}\right).
\end{align}
Taylor expanding $A_\mu$ in the first term to order $\delta$ gives
\begin{align}
\delta V = \hat{P}\exp \left( ig\int_0^1 \dd t \left( A_\mu(s) + \frac{ \D A_\mu}{\D s^\nu} \delta s^\nu\right)\left(\frac{ \dd s^\mu}{\dd t} + \frac{ \dd \delta s^\mu}{\dd t}\right)\right) - \hat{P} \exp\left(ig\int_0^1\dd t A_\mu(s) \frac{ \dd s^\mu}{\dd t}\right).
\end{align}
Multiplying out the parentheses and ignoring the term of order $\delta^2$, we have
\begin{align}
\delta V = \hat{P}\exp \left( ig\int_0^1 \dd t \left( A_\mu(s) \frac{ \dd s^\mu}{\dd t} + A_\mu(s) \frac{ \dd \delta s^\mu}{\dd t} + \frac{ \D A_\mu}{\D s^\nu} \delta s^\nu \frac{ \dd s^\mu}{\dd t} \right)\right) - \hat{P} \exp\left(ig\int_0^1\dd t A_\mu(s) \frac{ \dd s^\mu}{\dd t}\right).
\end{align}
Notice we can rewrite the $A_\mu(s) \frac{ \dd s^\mu}{\dd t}$ as $\frac{ \dd}{\dd t} \left( A_\mu(s) \delta s^\mu\right) - \frac{ \D A_\mu}{\D s^\nu} \frac{ \dd s^\nu}{\dd t} \delta s^\mu.$ [We used the chain rule to write $\frac{ \dd A_\mu}{\dd t} = \frac{ \D A_\mu}{\D s^\nu} \frac{ \dd s^\nu}{\dd t}$.] This is just integration by parts but keeping the boundary term.
\begin{align}
\delta V = \hat{P} \exp\left(ig \int_0^1 \dd t\left( A_\mu(s) \frac{ \dd s^\mu}{\dd t} + \frac{ \dd}{\dd t} (A_\mu(s) \delta s^\mu) + \left(\frac{ \D A_\mu}{\D s^\nu} - \frac{ \D A_\nu}{\D s^\mu}\right) \delta s^\nu \frac{ \dd s^\mu}{\dd t}\right)\right) - \hat{P} \exp\left(ig\int_0^1\dd t A_\mu(s) \frac{ \dd s^\mu}{\dd t}\right).
\end{align}
Now, notice that the path ordering operator will enfore ordering. Thus all matrices inside the exponent will commute. So we can insert commutators without changing the result, since any commutator is zero. Thus let's add an $ig[A_\mu,A_\nu]$ to complete the field strength tensor $F_{\nu\mu}$.
\begin{align}
\delta V = \hat{P} \exp\left(ig \int_0^1 \dd t\left( A_\mu(s) \frac{ \dd s^\mu}{\dd t} + \frac{ \dd}{\dd t} (A_\mu(s) \delta s^\mu) + F_{\mu\nu} \delta s^\mu \frac{ \dd s^\nu}{\dd t}\right)\right) - \hat{P} \exp\left(ig\int_0^1\dd t A_\mu(s) \frac{ \dd s^\mu}{\dd t}\right).
\end{align}
In order to simplify what follows, let's define $F(t) \equiv igA_\mu(s)\frac{ \dd s^\mu}{\dd t}$ and $G(t) \equiv ig\frac{ \dd}{\dd t} (A_\mu(s) \delta s^\mu) + igF_{\mu\nu} \delta s^\mu \frac{ \dd s^\nu}{\dd t}$. Then we have that $F(t)$ is $\mathcal{O}(\delta^0)$ and $G(t)$ is $\mathcal{O}(\delta)$.
\begin{align} \label{deltav}
\delta V = \hat{P} \exp\left( \int_0^1 \dd t\left( F(t) + G(t)\right)\right) - \hat{P} \exp\left(\int_0^1\dd t \,F(t)\right).
\end{align}
Expanding the first exponential, we then have
\begin{align}
\hat{P}\left( 1+ \int \dd t(F(t) + G(t)) + \frac12\int \dd t^\prime \dd t(F(t) + G(t))(F(t^\prime) + G(t^\prime)) + ...\right)
\end{align}
Then let's cancel any terms of order $\delta^2$ or higher. This means that we can only keep terms with $F(t)'s$ only and terms with only a single $G(t)$. Due to the path ordering operator, we can treat all the matrices inside the parentheses as commuting. Thus, it simplifies to:
\begin{align}
\hat{P}\left( \exp\left( \int \dd t F(t)\right) + \int \dd t G(t)\left( 1+ \int \dd t^\prime F(t^\prime) + \frac12\int \dd t^\prime \dd t^{\prime\prime} F(t^{\prime})F(t^{\prime\prime}) + ...\right)\right).
\end{align}
Putting this all back as the first term in Eq. \ref{deltav}, we see that we're left with:
\begin{align}
\delta V = \hat{P} \int \dd t\, G(t) \exp\left( \int \dd t^\prime\, F(t^\prime)\right).
\end{align}
Now we can put back in the definition of $F(t)$ and $G(t)$:
\begin{align}
\delta V = \hat{P} \int \dd t\, ig\left( \frac{ \dd}{\dd t} (A_\mu(s)\delta s^\mu) + F_{\mu \nu} \delta s^\mu \frac{ \dd s^\nu}{\dd t}\right)\exp\left(ig \int \dd t^\prime A_\sigma(s) \frac{ \dd s^\sigma}{\dd t^\prime}\right).
\end{align}
Notice that the exponential at the end here is just the Wilson line itself, except that now it needs to be path ordered along with the first integral over $t$. Let's now split the first integral into its two terms. The first integral is simply
\begin{align}
\int \dd t \frac{ \dd}{\dd t}(A_\mu(s)\delta s^\mu) = A_\mu(s(1))\delta s^\mu(1) - A_\mu(s(0))\delta s^\mu(0).
\end{align}
Which can be rewritten in another integral form:
\begin{align}
\int \dd t\, A_\mu(s)\delta s^\mu ( \delta(t-1) - \delta(t)).
\end{align}
We are being sloppy with the bounds on the integrations, but we could always go back to the beginning and define the integral to be over the range $t \in (-\infty,\infty)$ with appropriate step functions to fix the bounds. Now we have
\begin{align}
\delta V = \hat{P}ig \int \dd t\left(A_\mu(s) (\delta(t-1) - \delta(t)) + F_{\mu \nu} \frac{ \dd s^\nu}{\dd t}\right) \exp\left(ig\int \dd t^\prime \, A_\sigma(s)\frac{ \dd s^\sigma}{\dd t^\prime}\right) \delta s^\mu.
\end{align}
Comparing with Eq. \ref{funcderiv} (Sorry, I'm not sure how to get hyperref to work here), we can now read off the functional derivative.
\begin{align}
\frac{ \delta V}{\delta s^\mu (t)} = \hat{P}\,ig\left[ A_\mu(s)(\delta(t-1) - \delta(t)) + F_{\mu \nu} \frac{ \dd s^\nu}{\dd t}\right] \exp\left(ig\int \dd t^\prime \, A_\sigma(s) \frac{ \dd s^\sigma}{\dd t^\prime} \right).
\end{align}
We can now insert this result into Eq. \ref{Dform}(sorry again!) to get the partial derivative of the Wilson line:
\begin{align}
\frac{ \D V(z,y)}{\D x^\lambda} = \hat{P}\,ig\int \dd t \left[ A_\mu(s)(\delta(t-1) - \delta(t)) + F_{\mu \nu} \frac{ \dd s^\nu}{\dd t}\right] \exp\left(ig\int \dd t^\prime \, A_\sigma(s) \frac{ \dd s^\sigma}{\dd t^\prime} \right) \frac{ \D s^\mu}{\D x^\lambda}
\end{align}
Here we can make some simplifications. First, we can separate the terms and perform the $t$ integral in two of the terms using the delta functions.
\begin{align}
\frac{ \D V(z,y)}{\D x^\lambda} =& \hat{P}\,ig A_\mu(z)\frac{ \D z^\mu}{\D x^\lambda}\exp\left( \int \dd t^\prime A_\sigma(s) \frac{ \dd s^\sigma}{\dd t^\prime}\right) - \hat{P}\,ig A_\mu(y) \frac{ \D y^\mu}{\D x^\lambda} \exp\left( \int \dd t^\prime A_\sigma(s) \frac{ \dd s^\sigma}{\dd t^\prime}\right)\\ &+\hat{P}\,ig \int \dd t\, F_{\mu \nu} \frac{ \dd s^\nu}{\dd t} \frac{ \D s^\mu}{\D x^\lambda} \exp\left( ig\int \dd t^\prime A_\sigma(s) \frac{ \dd s^\sigma}{\dd t^\prime}\right).
\end{align}
Since $z$ is already the end point of the path, in the first term we can pass the path ordering on to the exponential, thus turning it into the Wilson line itself. In the second term, since $y$ is the starting point, we move $A_\mu(y)$ to the right of the exponential, and then put the path ordering operator on the exponential alone, thus also creating the Wilson line. In the third term, let's "un-parametrize" the path so that $s(t) \rightarrow s$ and $s(t^\prime) \rightarrow s^\prime$.
\begin{align}
\frac{ \D V(z,y)}{\D x^\lambda} =ig A_\mu(z) V(z,y) \frac{ \D z^\mu}{\D x^\lambda} -ig V(z,y)A_\mu(y)\frac{ \D y^\mu}{\D x^\lambda} + \hat{P}\,ig\int_y^z F_{\mu \nu}(s) \dd s^\nu \exp\left(ig \int_y^z A_\sigma(s^\prime) \dd s^{\sigma^\prime}\right)\frac{ \D s^\mu}{\D x^\lambda}.
\end{align}
In the third term, we can account for the path ordering by splitting the exponential into two pieces, the part of the path before $s$, and the part of the path after $s$. So,
\begin{align}
&\hat{P} \,ig\int_y^z F_{\mu \nu}(s) \dd s^\nu \exp\left(ig \int_y^z A_\sigma(s^\prime) \dd s^{\sigma^\prime}\right) \frac{ \D s^\mu}{\D x^\lambda} \\ &= ig\int_y^z \dd s^\nu \left[ \hat{P}\exp\left(ig\int_s^z A_\sigma(s^\prime) \dd s^{\sigma^\prime}\right)\right] F_{\mu \nu}(s) \frac{ \D s^\mu}{\D x^\lambda} \left[ \hat{P}\exp\left(ig\int_y^s A_\tau(s^\prime)\dd s^{\tau^\prime}\right) \right]\\
&= ig\int_y^z \dd s^\nu\, V(z,s)F_{\mu \nu}(s) \frac{ \D s^\mu}{\D x^\lambda} V(s,y).
\end{align}
We can also write the first two terms in a simplified way:
\begin{align}
igA_\mu(z) V(z,y) \frac{ \D z^\mu}{\D x^\lambda} - igV(z,y)A_\mu(y)\frac{ \D y^\mu}{\D x^\lambda} = igV(z,s)A_\mu(s) \frac{ \D s^\mu}{\D x^\lambda} V(s,y)|^{s=z}_{s=y}.
\end{align}
Where we've exploited the fact that $V(z,z) = V(y,y) = 1$. Thus our final result is as follows.
\begin{align}
\frac{ \D V(z,y)}{\D x^\lambda} = igV(z,s)A_\mu(s) \frac{ \D s^\mu}{\D x^\lambda} V(s,y)|^{s=z}_{s=y} + ig\int_y^z \dd s^\nu\, V(z,s)F_{\mu \nu}(s) \frac{ \D s^\mu}{\D x^\lambda} V(s,y).
\end{align}
Where the first term is used if $x$ is either the start point or the end point for the Wilson line, and the last term is used if $x$ is an interior point of the Wilson line.
