Let's simplify slightly. Suppose the overdetermined $m\times n$ system $A$ admits the factorization $Q R = A$ with $R$ full rank. Let
$$A^\dagger = R^{-1} Q^T$$
be the linsolve
solution operator and observe that
$$I = R^{-1} Q^T Q R = A^\dagger A.$$
Differentiating the above, we have
$$ 0 = \frac{\partial (A^\dagger A)}{\partial g} = \frac{\partial A^\dagger}{\partial g} A + A^\dagger \frac{\partial A}{\partial g} . $$
We cannot solve for $\frac{\partial A^\dagger}{\partial g}$ algebraically because it is not uniquely determined in this form, but we can compute
$$ \frac{\partial A^\dagger}{\partial g} Q Q^T = - A^\dagger \frac{\partial A}{\partial g} R^{-1} Q^T = - R^{-1} Q^T \frac{\partial A}{\partial g} R^{-1} Q^T = - A^\dagger \frac{\partial A}{\partial g} A^\dagger $$
where $QQ^T$ is the identity on the range of $A$.
Edit
If we need the rest, $\frac{\partial A^\dagger}{\partial g} (I - Q Q^T)$, we can formally differentiate the QR factorization,
$$ \frac{\partial A^\dagger}{\partial g} = \frac{\partial R^{-1}}{\partial g} Q^T + R^{-1} \frac{\partial Q^T}{\partial g},$$
and decompose into the parts in the range of $A$ (which yields expressions equivalent to that derived above) and in the left null space of $A$,
$$\begin{align} \frac{\partial A^\dagger}{\partial g} (I - Q Q^T) &= \left[ \frac{\partial R^{-1}}{\partial g} Q^T + R^{-1} \frac{\partial Q^T}{\partial g} \right] (I - Q Q^T) \\
&= R^{-1} \frac{\partial Q^T}{\partial g} (I - Q Q^T) .
\end{align} $$
Note that only perturbations $dg$ that change $Q$ can contribute to this second term. We can rewrite in terms of $\frac{\partial A}{\partial g}$ by
$$\begin{align}
\frac{\partial A^\dagger}{\partial g} (I - Q Q^T) &= R^{-1} R^{-T} \left[ R^T \frac{\partial Q^T}{\partial g} + \frac{\partial R^T}{\partial g} Q^T \right] (I - Q Q^T) \\
&= R^{-1} R^{-T} \frac{\partial A^T}{\partial g} (I - Q Q^T) \\
&= (A^T A)^{-1} \frac{\partial A^T}{\partial g} (I - A A^\dagger)
\end{align}$$
where the last line is the second term in Stefano's derivation, though the QR form
$$\frac{\partial A^\dagger}{\partial g} = R^{-1} \left[ -Q^T \frac{\partial A}{\partial g} R^{-1} Q^T + R^{-T} \frac{\partial A^T}{\partial g} (I - Q Q^T) \right]$$
is preferable to compute with for numerical stability and because it reuses more computation and operates in a smaller space.
You can apply similar arguments to rank-deficient and under-determined cases.