We analyze strongly interacting Fermi gases in the unitary regime by considering the generalization to an arbitrary number N of spin-1/2 fermion flavors with Sp(2N) symmetry. For N=infty this problem is exactly solved by the BCS-BEC mean-field theory, with corrections small in the parameter 1/N. The large-N expansion provides a systematic way to determine corrections to mean-field predictions, allowing the calculation of a variety of thermodynamic quantities at (and in the proximity to) unitarity, including the energy, the pairing gap, and upper-critical polarization (in the case of a polarized gas) for the normal to superfluid instability. For the physical case of N=1, among other quantities, we predict in the unitarity regime, the energy of the gas to be xi=0.28 times that for the non-interacting gas and the pairing gap to be 0.52 times the Fermi energy.