We introduce a new fermionic variational wavefunction, generalizing the Bardeen-CooperSchrieffer (BCS) wavefunction, which is suitable for interacting multi-species spinful systems and sustaining superfluidity. Applications range from quark matter to the high temperature superconductors. A wide class of Hamiltonians, comprising interactions and hybridization of arbitrary momentum dependence between different fermion species, can be treated in a comprehensive manner. This is the case, as both the intra-species and the inter-species interactions are treated on equally rigorous footing, which is accomplished via the introduction of a new quantum index attached to the fermions. The index is consistent with known fermionic physics, and allows for heretofore unaccounted fermion-fermion correlations. We have derived the finite temperature version of the theory, thus obtaining the renormalized quasiparticle dispersion relations, and we discuss the appearance of charge and spin density wave order.We present numerical solutions for two electron species in 2 dimensions. Based on these solutions, we show that, for equivalent spin up and down fermions, the Fermi occupation factor (per spin) equals 1/2 deep in the Fermi sea. This constitutes a unique experimental prediction of the theory, both for the normal and superfluid states. Interestingly, this result, obtained in the thermodynamic limit, is consistent with Fermi occupation factor (in-)equalities for finite systems of electrons, derived (in a different context) by Borland and Dennis, J. Phys. B 5, 7 (1972) and by Altunbulak and Klyachko, Commun. Math. Phys. 282, 287 (2008).