Assuming you refer to the ZOH-discretization as shown in this figure from the mathworks site
the relationship between the continuous-time signals and the discrete-time signals can be derived as follows. The signal $u(t)$ is given by
$$u(t)=\sum_ku[k]g(t-kT)\tag{1}$$
where $g(t)$ is a rectangular impulse response (constant in the interval $t\in [0,T]$, zero anywhere else) representing the ZOH, and $T$ is the sampling period. Let $h(t)$ be the impulse response corresponding to the transfer function $H(s)$. Furthermore, let $f(t)$ be the convolution of $g(t)$ and $h(t)$:
$$f(t)=(g\star h)(t)\tag{2}$$
Then the continuous-time output signal $y(t)$ is given by
$$\begin{align}y(t)&=(u\star h)(t)\\
&=\int_{-\infty}^{\infty}u(\tau)h(t-\tau)d\tau\\&=\sum_ku[k]\int_{-\infty}^{\infty}g(\tau-kT)h(t-\tau)d\tau\\&=\sum_ku[k]f(t-kT)\tag{3}\end{align}$$
where I've used Equations $(1)$ and $(2)$. From $(3)$ the discrete-time output signal is easily obtained as
$$y[n]=y(nT)=\sum_ku[k]f((n-k)T)=\sum_ku[k]f[n-k]\tag{4}$$
with $f[n]=f(nT)$. So $y[n]$ is simply the discrete-time convolution of the input signal $u[n]$ and the total impulse response $f[n]$.
From this result it immediately follows that linearity must be satisfied, simply because convolution is a linear operation:
$$f(t)=\alpha_1(g\star h_1)(t)+\alpha_2(g\star h_2)(t)=(g\star (\alpha_1h_1+\alpha_2h_2))(t)\tag{5}$$
However, the (ZOH-)discretization of the concatenation of two transfer functions is not equivalent to the concatenation of the two discretized transfer functions. Referring to the figure above, the first case is equivalent to replacing $H(s)$ by the concatenation of $H_1(s)$ and $H_2(s)$. The other case involves concatenating two complete systems as shown in the figure, the first with $H(s)=H_1(s)$, the second with $H(s)=H_2(s)$. The difference between the two cases is that in the first case you only have one ZOH, whereas in the second one you get two ZOHs. The equivalent discrete-time impulse responses are
$$f[n]=(g\star h_1\star h_2)(nT)\tag{6}$$
and
$$f[n]=(g\star h_1)\star (g\star h_2)(nT)\tag{7}$$
which are generally not identical. The first case in Eq. $(6)$ is a ZOH discretization of the total transfer function $H(s)=H_1(s)H_2(s)$, whereas the second case in Eq. $(7)$ corresponds to a first-order hold (FOH) discretization of the total transfer function.
cont2discrete
withmethod='zoh'
on several combinations of LTI systems, and it actually seems to hold exactly in every case that I try. The discretization may be nonlinear in terms of the state-space matrices, but so is the mapping from state-space to t.f. ($F(s) = C(sI - A)^{-1}B + D$) -- and so the two nonlinearities may 'cancel'. However I can't find a proof of this claim anywhere online. It may just follow from the exactness of the discretization, and the fact these systems are linear. $\endgroup$