The present paper deals with the problem of constructing discrete counterparts of conjugate harmonic functions in the scope of hypercomplex analysis. Specifically, let us denote by \(\Delta _h\) the discrete Laplacian, by \(D_h\) the discrete Dirac operator such that \((D_h)^2=-\Delta _h\) , and by \(\partial _s\) a time-derivative. Our approach starts with the correspondence between the null solutions \(\varvec{w}=\varvec{u}-\textbf{e}_0\varvec{v}\) of the Dirac type operator \(\textbf{e}_0\partial _s+D_h\) and the pair of solutions \((\varvec{u},\varvec{v})\) of the following semidiscrete Cauchy-Riemann type system: \(\begin{aligned} {\left\{ \begin{array}{ll} \partial _s\varvec{u}({\varvec{x}},s)=-D_h \varvec{v}({\varvec{x}},s) & \\ \partial _s\varvec{v}({\varvec{x}},s)=-D_h \varvec{u}({\varvec{x}},s)& \end{array}\right. },&({\varvec{x}},s)\in h{{\mathbb {Z}}}^n \times (0,\infty ). \end{aligned}\) Hereby, \(\textbf{e}_0\) stands for a Clifford basis element, satisfying \(\textbf{e}_0^2=-1\) and \(\textbf{e}_0D_h+D_h\textbf{e}_0=0\) . In a nutshell, we prove that the solutions of such system can be characterized in terms of the null solutions of the semidiscrete Laplacian \(\partial _s^2+\Delta _h\) . Afterwards, we show that a similar formulation on discrete space-time lattices arises from the aforementioned semidiscrete formulation on \(h{{\mathbb {Z}}}^n \times (0,\infty )\) by means of sampling using Bessel functions of the first kind. These results provide us the building blocks to characterize the discrete analogue of the Riesz-Hilbert transform in terms of the ’boundary behavior’ of discrete conjugate harmonic functions, generated from operator semigroups of Poisson type. We end this paper by discussing some open problems of future research.