# A closed-form formula for Cov(X,Y) when X and Y are normal random vectors?

I cannot figure out one step given in my textbook, [1] Mixed Models. Theory and Applications by E. Demidenko.

I study the Linear Mixed Effects (LME) model in the following form: $$\mathbf{y_i}=\mathbf{X_i\beta}+\mathbf{Z_ib_i} + \varepsilon _i.$$ It is assumed that $$\mathbf{b_i} \sim \mathit{N}(\mathbf{0},\sigma ^2\mathbf{D})$$ and, eventually, that the LME model with normally distributed random variables can be written in marginal form as $$\mathbf{y_i} \sim \mathit{N}(\mathbf{X_i\beta},\sigma ^2\mathbf{(I+Z_iDZ'_i)}) \tag{2.13}$$

What puzzles me is the part about $\operatorname{cov}(\mathbf{b_i},\mathbf{y_i})$ ([1], p. 148.):

Assuming normal distribution, from model (2.13) it follows that $\operatorname{cov}(\mathbf{b_i},\mathbf{y_i}) = \sigma ^2\mathbf{DZ'_i}$.

I guess that this follows some known fact about covariance between two normal random vectors (?) but I cannot figure out where from we do have this $\operatorname{cov}(\mathbf{b_i},\mathbf{y_i})$ formula.