Consider an elliptic parameter $k$; we introduce a family of $Z^u$-Dirac operators $(\mathsf{K}(u))_{u\in\Re(\mathbb{T}(k))}$, relate them to the $Z$-massive Laplacian of [BdTR17b], and extend to the full $Z$-invariant case the results of Kenyon [Ken02] on discrete holomorphic and harmonic functions, which correspond to the case $k=0$. We prove, in a direct statistical mechanics way, how and why the $Z^u$-Dirac and $Z$-massive Laplacian operators appear in the $Z$-invariant Ising model, considering the case of infinite and finite isoradial graphs. More precisely, consider the dimer model on the Fisher graph ${\mathsf{G}}^{\scriptscriptstyle{\mathrm{F}}}$ arising from a $Z$-invariant Ising model. We express coefficients of the inverse Fisher Kasteleyn operator as a function of the inverse $Z^u$-Dirac operator and also as a function of the $Z$-massive Green function; in particular this proves a (massive) random walk representation of important observables of the Ising model. We prove that the squared partition function of the Ising model is equal, up to a constant, to the determinant of the $Z$-massive Laplacian operator with specific boundary conditions, the latter being the partition function of rooted spanning forests. To show these results, we relate the inverse Fisher Kasteleyn operator and that of the dimer model on the bipartite graph ${\mathsf{G}}^{\scriptscriptstyle{\mathrm{Q}}}$ arising from the XOR-Ising model, and we prove matrix identities between the Kasteleyn matrix of ${\mathsf{G}}^{\scriptscriptstyle{\mathrm{Q}}}$ and the $Z^u$-Dirac operator, that allow to reach inverse matrices as well as determinants.