Use of Itô isometry for correlation calculation

Question

When calculating the covariance of the Ornstein-Uhlenbeck process, the Wikipedia article applies implicitly the Itô isometry with the fact of non-overlapping independent increments of the Wiener process in the following step: $$ \sigma^2 e^{-\theta (s+t)}\,\mathbb{E} \left[ \int_0^s e^{\theta u}\, \mathrm{d}W_u \int_0^t e^{\theta v}\, \mathrm{d}W_v \right] = \frac{\sigma^2}{2\theta} \, e^{-\theta (s+t)}(e^{2\theta \min(s,t)}-1) $$ Could someone explain me how is this exactly done, step by step? More generally, how does one apply the Itô isometry to these two processes? $$ \mathbb{E} \left[ \int_0^s f(u)\,\mathrm{d} W_u \int_0^t g(v) \,\mathrm{d} W_v \right] $$