You're correct, in an undergrad statistical mechanics course we are taught that we can study metals as a free electron gas, and indeed it gives you (roughly) the correct physical quantities, like specific heat and compressibility. There are two problems with this pictures. First, the electrons are contained in a lattice. Bloch's theorem tells us that free electrons in a lattice behave as if they we're in homogeneous space, except for some amplitudes, so we're fine. Therefore let's ignore the lattice (and the phonons too) for simplicity, the so called jellium model. We now have to explain the second problem, namely that the electron-electron coulomb interaction will not change radically the free electron conclusions. This is what Fermi's liquid theory is for, as you've noticed.
Now, all equilibrium thermodynamical properties are dictated by the ground state (or vacuum in QFT language) in zero temperature, or by the ground state plus excited states with statistical weight in non-zero temperature. Let me focus first on the $T=0$ regime and then we can discuss the more general case. In the free electron model the ground state is just the Fermi surface, and the excited states are electron-hole pairs created in the surface. What about the interacting case? The key here is Gell-Mann-Low Theorem, which says (roughly) that if the interaction is weak we can connect the ground state of the non-interacting theory adiabatically with the interacting one. See, in the interacting case the coulomb repulsion will force some electrons out of the Fermi surface in a very complicated way, making the new ground state difficult to analyse with usual perturbation theory. The excitations are indeed excited states full stop, is just that without the Fermi surface as reference we don't really know who they are. The theorem allows us to calculate everything in terms of the non-interacting system using the perturbative calculations for the Green's functions. Why care about excitations at all? Well, according to linear response theory, quantities such as magnetic susceptibility and conductivity can be framed as a scattering process between excited states. Using Green's functions we can calculate everything based on the non-interacting system which is great. So it's important to look at this thing in order to have measurable physical quantities.
Now, in the free case the Green's function (in momentum-frequency space) is
$G(\omega,k)=\frac{1}{\omega-\epsilon_k}$
where $\epsilon_k$ is the energy of the state with momentum $k$, and I'm ignoring the chemical potential for simplicity. The spectral function $A(\omega,k)=-\pi^{-1}\mathcal{Im}\{G\}$ will be just a Dirac delta, which indicates a particle behavior, that is energy with a definite momentum. Now, the interacting Green's function will be
$G(\omega,k)=\frac{Z_k}{\omega-E_k+i\gamma_k}$+ some background noise
Where $E_k$ is the energy evaluated in perturbation theory. The background noise will appear as a non-zero signal in every momentum, reflecting the lots of different electron-hole pairs, but the first part will give a spectral function which looks like a gaussian package centered at $E_k$ with a width $\gamma_k$. This looks almost like the particle distribution, hence quasi-particle, except that the width $\gamma_k$ implies that this quasi-particles eventually decay, by interaction, into the background noise. But if the width is very small then the decay times can be very long, and almost everything will seem like the free electron gas. Note that $E_k=v_F(k-k_F)$ near the Fermi surface, but $\gamma_k\propto (k-k_F)^2$, where $k_F$ is the Fermi momentum, and $v_F$ the Fermi velocity, so that for low excitations it is true that the width vanishes faster than the energy, therefore justifying the free electron gas. The weight $Z_k<1$ will tell you how much of the system behaves in this way.
It should be intuitive that this quasi-particles are still fermionic, but in any case in a finite temperature formalism, such as Matsubara's, we can strictly prove that they obey Fermi-Dirac statistics, and finite temperature proceeds smoothly from there.
So the ingenious part is that, if the perturbation is small, one can use the whole machinery from QFT to evaluate everything in terms of the non-interacting case, and the computations show that in the general case you get Green's functions that almost look like free particle ones, except for some renormalized energies and a decay probability.
Just for completeness, note that for high momentum excitations the width grows, so somewhere all those approximations break. Also, it could be the case that the coulomb repulsion is not a small perturbation, invalidating Gell-Mann-Low Theorem. In the strong interacting regime some different thing might appear, such as a Mott insulator.
