It might help to consider a scalar example to start. If you are familiar with statistical mechanics, consider the partition function of a system that consists of three spin 1/2 particles (so there are 8 states total), and let us suppose that the Hamiltonian can be written
\begin{align*}
H &= -h\sum_{i=0}^2\sigma_i-J\sum_{i=0}^2\sigma_i\sigma_{\overline{i+1}},\quad \overline{i}\equiv i\mod 3\\
&\equiv H_0(\sigma_0,\sigma_1,\sigma_2) +H_{int}(\sigma_0,\sigma_1,\sigma_2)
\end{align*}
The partition function, $Z(h,J) = \sum_{\sigma_0,\sigma_1,\sigma_2}e^{- H(\sigma_0,\sigma_1,\sigma_2)}$ can be computed directly:
\begin{align*}
Z(h,J) = e^{3(h + J)} + 3e^{(h-J)} + 3e^{-(h-J)}+e^{-3(h+J)}
\end{align*}
Recall that in quantum mechanics, the derivative $-i\partial_x$ (often used to represent the momentum operator) can be viewed as a generator of translations: $\exp(-m\partial_x)\Psi(x) = \Psi(x-m)$. (This is essentially because $f(x+h)\approx f(x)+h\nabla f(x)$.) The exact same principles may be applied in the example above: to evaluate $Z(h,J-m)$, simply apply the operator $\exp(-m\partial_J)$ to $Z(h,J)$.
Now let's explore various ways of interpreting $\exp(-m\partial_J)$, or $\partial_J$ even more simply. Just to give a small taste the full set of available approaches (especially in lattice models), we'll start with a somewhat tangential example.
First, notice in
\begin{align*}
\partial_JZ(h,J) &= -\sum_\sigma M_{int}(\sigma)e^{-hM_0(\sigma) - JM_{int}(\sigma)}\\
&=\sum_\sigma (\sigma_0\sigma_1 + \sigma_1\sigma_2 + \sigma_2\sigma_0)e^{-h\sum_{i=0}^2\sigma_i-J\sum_{i=0}^2\sigma_i\sigma_{\overline{i+1}}}
\end{align*}
that the expression $\sigma_0\sigma_1+\sigma_1\sigma_2+\sigma_2\sigma_0$ would match $\frac{1}{2}(\sigma_0 + \sigma_1 + \sigma_2)^2$ if it weren't for the extra $\frac{1}{2}(\sigma_0^2+\sigma_1^2+\sigma_2^2)$ term. In fact, because the $\sigma_i$ are limited to $\{-1, 1\}$, the last expression ($\frac{1}{2}\sum_i\sigma_i^2$) is just a constant; $\frac{3}{2}$ in this case. More generally, however, we could introduce a new term to the Hamiltonian, $-P\sum_{i}\sigma_i^2\equiv P M_{sq}(\sigma)$, so that $M_{int}(\sigma) = \frac{1}{2}(M_0(\sigma)^2 - M_{sq}(\sigma))$. Hence,
\begin{align*}
\partial_J Z(h,J) = -\sum_\sigma \frac{1}{2}(M_0(\sigma)^2 - M_{sq}(\sigma))e^{-hM_0(\sigma)-JM_{int}(\sigma)-PM_{sq}(\sigma)}
\end{align*}
(where the right hand side is understood as being evaluated at $P = 0$.) Notice that $-M_{sq}(\sigma)e^{-PM_{sq}(\sigma)}$ can be replaced with $\partial_P e^{-PM_{sq}(\sigma)}$, and $M_0(\sigma)^2 e^{-hM_0(\sigma)}$ can be replaced with $\partial_h^2 e^{-hM_0(\sigma)}$. We have just derived the relation
\begin{align*}
\partial_J Z(h,J) = -\frac{1}{2}\partial_h^2 Z(h,J,P) - \partial_P Z(h,J,P)\Big|_{P=0}
\end{align*}
More possibilities present themselves if we allow the parameters $h$, $J$, $P$ etc. to assume different values for each spin (for example, the spins might be impurities positioned variously in a solid, and form an irregular triangle so that pairwise neighbor couplings are distinct.) In this case, the Hamiltonian would read
\begin{align*}
H(\sigma) = -\sum_{i=0}^2(h_i \sigma_i - J_i \sigma_i\sigma_{\overline{i+1}})
\end{align*}
with the partition function modified appropriately.
Instead of having to introduce an auxiliary $M_{sq}(\sigma)$ term, we can now simply write
\begin{align*}
\sigma_0\sigma_1 e^{-\sum_i h_i\sigma_i} = \partial_{h_0}\partial_{h_1} e^{-\sum_i h_i\sigma_i}
\end{align*}
and
\begin{align*}
\partial_J Z(h,J) = \sum_{\sigma} \sum_{i=0}^2 \partial_{h_i} \partial_{h_{\overline{i+1}}} e^{-H_0(\sigma, h) - H_{int}(\sigma, J)}.
\end{align*}
(where $h$ is now understood as a three component vector or tuple.)
You may recognize $\sum_{i=0}^2\partial_{h_i}\partial_{h_{\overline{i+1}}}$ as $M_{int}(\partial_h)$, and moreover that $M_{int}(\partial_h)$ is independent of the spin degrees of freedom and can be moved outside the sum: written explicitly,
\begin{align*}
\partial_J Z(h_0, h_1, h_2, J) = -M_{int}(\partial_{h_0}, \partial_{h_1}, \partial_{h_2}) Z(h_0, h_1, h_2, J).
\end{align*}
Finally, having established that the linear differential operators $\partial_J$ and $M_{int}(\partial_h)$ are equivalent, we can substitute one for the other in the exponential:
\begin{align*}
\exp(-m\partial_J) = \exp(-m M_{int}(\partial_h)),
\end{align*}
and hence,
\begin{align*}
Z(h,J-m) = \exp(-m M_{int}(\partial_h)) Z(h,J).
\end{align*}
