In Ref. Gaudin (Nuclear Physics 15, 89 1960), Gaudin established a generalized version of Wick’s theorem at finite temperature for systems described by a diagonal quadratic Hamiltonian. In this work, we extend Gaudin’s proof to statistical density operators generated by a general quadratic Hamiltonian, prior to diagonalization, as commonly encountered in interacting Bose systems. We show explicitly that the resulting thermal contractions depend on the Bogoliubov transformation parameters and on the quasiparticle spectrum, providing a transparent framework for evaluating thermal averages of arbitrary operator products. As an illustration of the usefulness of this generalized formulation, we evaluate particle-number fluctuations in a homogeneous, weakly interacting Bose gas below the Bose–Einstein condensation temperature. We show that particle-number fluctuations are generically non-Poissonian as a consequence of interaction-induced correlations and quantum depletion, with the Poissonian limit recovered only for an ideal Bose gas at zero temperature. Our results clarify the role of quadratic Hamiltonians and Bogoliubov quasiparticles in shaping fluctuation statistics and highlight the generalized Wick’s theorem as a controlled tool for computing correlation functions in interacting quantum gases.