We have computed pKa shifts for carboxylic residues of the serine protease inhibitor turkey ovomucoid third domain (residues Asp7, Glu10, Glu19, Asp27, and Glu43). Both polarizable and nonpolarizable empirical force fields were employed. Hydration was represented by the surface generalized Born and Poisson-Boltzmann continuum model. The calculations were carried out in the most physically straightforward fashion, by directly comparing energies of the protonated and deprotonated protein forms, without any additional parameter fitting or adjustment. Our studies have demonstrated that (i) the Poisson-Boltzmann solvation model is more than adequate in reproducing pKa shifts, most likely due to its intrinsically many-body formalism; (ii) explicit treatment of electrostatic polarization included in our polarizable force field (PFF) calculations appears to be crucial in reproducing the acidity constant shifts. The average error of the PFF results was found to be as low as 0.58 pKa units, with the best fixed-charges average deviation being 3.28 units. Therefore, the pKa shifts phenomena and the governing electrostatics are clearly many-body controlled in their intrinsic nature; (iii) our results confirm previously reported conclusions that pKa shifts for protein residues are controlled by the immediate environment of the residues in question, as opposed to long-range interactions in proteins. We are confident that our confirmation of the importance of explicit inclusion of polarization in empirical force fields for protein studies will be useful far beyond the immediate goal of accurate calculation of acidity constants.