How to derive covariance between Y and random effect in linear mixed model In Searle et al.'s book Variance Components, the covariance between the response variable Y and the random effect u is described as follows (highlighted with red rectangle):

I am not sure why cov(y, u^T) is ZD. I want to see proof of this using expectation calculations beginning with something like cov(y, u^T) = E(y*u^T) - E(y)E(u^T)... or cov(y, u^T) = E[(y - E(y)) (u^T - E(u^T))]. I just do not know how to proceed. Thank you.

Fllowing Frank's answer, I also wish to calculate cov(u,y)
cov(u, y) = cov(u, Xβ+Zu+e)
= cov(u, Zu)
= E[uuTZT] - E[u]E[uTZT] =E[uuT]ZT-E[u]E[uT]ZT OR =E[uuT]ZT-E[u]ZTE[uT]
Following the equation before OR, I could proceed and get DZT. However, why is the equation after OR incorrect?

First, since covariance is bilinear and since $X\beta$ is constant and $e$ is independent noise, i.e. independent of $u$, we have:
$$
\begin{align}
\operatorname{cov}(y, u) &= \operatorname{cov}(X\beta + Zu + e, u)\\
          &= \operatorname{cov}(X\beta, u)+\operatorname{cov}(Zu, u) + \operatorname{cov}(e, u)\\
          &= 0 + \operatorname{cov}(Zu, u) + 0\\
          &= \operatorname{cov}(Zu, u).
\end{align}
$$
Now, by the definition of $\operatorname{cov}$ given by you and, because the expectation is linear, it follows:
$$
\begin{align}
\operatorname{cov}(y, u) &= \operatorname{cov}(Zu, u)\\
    &= E[Zuu^T] - E[Zu]E[u^T]\\
    &= ZE[uu^T] - ZE[u]E[u^T]\\
    &= Z(E[uu^T] - E[u]E[u^T])\\
    &= Z \operatorname{cov}(u)\\
    &= Z D.\\
\end{align}
$$
