If $|psirangle$ is a free fermionic state, why does its reduced density matrix $text{Tr}_C(|psirangle langle psi|)$ also obey Wick's theorem?

I have recently been trying to understand this paper. So far I understand why, given a free fermionic state $|psirangle$, it is fully characterised by its 2-point correlation matrix (i.e. obeys Wicks’s theorem). I also understand why if its reduced density matrix $rho = text{Tr}_C(|psiranglelangle psi|)$ obeys Wick’s theorem, that is $C_{ij} = text{Tr}(rho c^{dagger}_ic_j)$ fully characterises $rho$, then $rho$ must be gaussian (or the exponential of a free fermionic Hamiltonian, i.e. $rho sim exp(-sum_{ij}h_{ij}c^{dagger}_ic_j)$ ).

But why does $rho$ in the first place need to obey Wick’s theorem?

I was thinking about it and all I could come up with was that $rho sim sum_n a_n |tildepsi_nrangle langle tildepsi_n|$ where $|tilde psi_nrangle$ is the part of eigenstate $|psi_nrangle$ not in subsystem $C$. However, from my previous questions and specially with the help of @NorbertSchuch, I have gathered that in general that form for $rho$ is not gaussian $iff$ doesn’t obey Wick’s theorem. So I’m certainly missing something!