The answer by DaftWullie is useful for more general Hamiltonians, but since you are working specifically with the Ising Hamiltonian, then you do not need auxiliary qubits or perturbations which are two of the worst disadvantages of using perturbative gadgets. If you use perturbative gadgets, you might require millions of extra qubits, and coupling strengths between those qubits spanning several orders of magnitude.
Since the Ising Hamiltonian is diagonal (only being comprised of Pauli-$Z$ operatiors), you can simply use classical quadratization techniques, which have been used in the operations research community for decades.
Probably the simplest way (but certainly not the most sophisticated or efficient way) to reduce the Hamiltonian $Z_1Z_2Z_3Z_4$ is to first convert the $Z$ operators (which have eigenvalues +1 and -1) to $B$ operators (which describe precisely the same physics but with the energies of all states shifted up by a constant, or the "zero" of energy moved by up by one unit):
\begin{align}
\tag{1}
B &\equiv \frac{1}{2}\left(Z + I\right) \\
&= \begin{pmatrix} 1 & 0 \\ 0 & 0 \end{pmatrix}.\tag{2}
\end{align}
Then the ground state of a Hamiltonian such as the following (with $\alpha,\beta$ and $\gamma$ being whatever coupling-constants you wish):
$$
H = \alpha B_1B_2B_5 + \beta B_1B_6B_7B_8 - \gamma B_3B_4\tag{3},
$$
can be found by minimizing the following three 2-local Hamiltonians:
\begin{align}
H_{0} &= -\gamma B_3B_4\tag{4}\\
H_{10} &= \alpha B_2B_5 - \gamma B_3 B_4 \tag{5}\\
H_{11} &= \alpha B_2B_5 + \beta B_7 B_8- \gamma B_3 B_4. \tag{6}
\end{align}
$H_0$ comes from setting $B_1 = 0$, while $H_{10}$ and $H_{11}$ come from setting $\{B_1,B_6\} = \{I,0\}$ and $\{I,I\}$ respectively. If you minimize each of the three Hamiltonians, then check to see which configuration gives the overall minimum energy, you will get precisely the ground state of the original 4-local Hamiltonian by only minimizing three (much simpler) 2-local Hamiltonians.
In the literature this has been called "split-reduction", and notice that the cost of doing this was sub-exponential in the number of splits: in this case we need to minimize three Hamiltonians after doing 3 splits, instead of $2^2=4$.
Alternatively if you used the perturbative gadget in DaftWullie's answer, you would end up with an incredibly more complicated Hamiltonian that not only now has $X$ operators (turning it from stoquastic to non-stoquastic, which makes it significantly harder to minimize whether on quantum annealers or classical simulators), but you'd also need 9 extra qubits: Notice that $Z_1Z_2Z_3$ became a Hamiltonian involving $\{Z_1,Z_2,Z_3,Z_A,Z_B,Z_C\}$ which means 3 extra qubits to quadratize the cubic term, which would be needed twice for the $B_1B_6B_7B_8$ term and again once for the $B_1B_2B_5$ term, leading to 9 extra qubits in total. Increasing the size of the Hilbert space by a factor of $2^9$ is much worse than minimizing 3 far simpler Hamiltonians that are in smaller Hilbert spaces than the original Hamiltonian.
