A comparative study of tetraheme cytochrome c3 molecules from several species was carried out using recently developed theoretical methods based on continuum electrostatics. The binding joint equilibrium of electrons and protons was simulated, revealing the complete thermodynamic aspects of electron-proton coupling in these molecules. The method yields excellent accuracy in terms of midpoint potentials, giving the correct reduction orders in all molecules examined, except for one heme site. The coupling between electrons and protons is shown to be present and significant at physiological pH in all cases. This phenomenon, known as the redox-Bohr effect, though of thermodynamic nature, is shown to have an intrinsic "dynamic" character at the molecular level (in the sense of the empty/occupied fluctuations at the microscopic level), with the binding states of redox and protonatable sites displaying both correlated averages and correlated fluctuations. The protonatable sites more directly involved in the redox-Bohr effect are identified using, among other properties, the statistical correlation between pairs of sites, which automatically reflects indirect effects mediated by other sites. Several sites are identified in this analysis. Propionate D of heme I seems to be the most interesting, generally showing a high correlation not only with its own heme, but also with heme II, corresponding to an indirect stabilization of the reduced forms of both hemes. Other interesting sites are the free histidines of two of the cytochromes and propionate D of heme IV, the latter being potentially associated with redox-induced structural changes. Among the set of cytochromes c3 analyzed in this study, significant differences are observed for several properties of the acidic cytochrome included in the set, from Desulfovibrioafricanus, supporting the hypothesis of a different functional role.