According to Wikipedia, the transitional probability of the Ornstein–Uhlenbeck process
$$ \mathrm dx_t = -\theta(\mu-x_t)\,\mathrm dt + \sigma\,\mathrm dW_t $$
is
$$ P(x, t\,|\,x_0, t_0) = \sqrt{-\frac{\theta}{\sigma^2(1-e^{-2\theta(t-t_0)})}}\exp\left(-\frac{\theta}{\sigma^2}\frac{(x-x_0e^{-\theta(t-t_0)})^2}{\sigma^2(1-e^{-2\theta(t-t_0)})}\right) $$
i.e.
$$ x_t\,|\,x_0 = N\left(x_0e^{-\theta(t-t_0)}, -\frac{1}{\theta}\sigma^2(1-e^{-2\theta(t-t_0)})\right) $$
without any proof. It mentions that this is derived from the Fokker–Planck representation, but would there be a way to derive this directly? Perhaps from the formal solution?
$$ x_t - x_0 e^{-\theta(t-t_0)} = \mu(1-e^{-\theta(t-t_0)}) + \sigma\int_{t_0}^t e^{-\theta(t-s)}\,\mathrm dW_s $$
Any references to literature with details would also be appreciated.