The rapid accumulation of entries in protein databases makes it possible to carry out large-scale sequence analysis and to make discoveries that cannot be sufficiently supported with a small dataset. Multiple sequence alignment is a representative method for identifying positional variations in related protein sequences by working with sizable entries. Despite the convenience of operation brought about by the straightforward comparison, the widely employed multiple sequence alignment seldom considers the covariation of two or more sites in a protein sequence. Such a simplification inevitably sacrifices the chance to discover knowledge that can better reflect heritable variations in the long history of evolution. The statistical coevolution analysis method goes further by focusing on the mutation rate of two positions to identify the functional correlation, based on which the independent components (IC) were then clustered into independent sections (sector). Using this method, we analyzed the metal-binding specificity of superoxide dismutases (SOD) from different families. The results showed that the residues coordinating metal ions could be clustered into a sector. However, the SODs of different families exhibited different levels of correlation, reflected by the order of the sectors: sector 1 for SODs from families PF09055 (NiSOD, based on the Pfam database) and PF00081(MnSODs or FeSODs), and sector 2 for PF00080 (Cu/ZnSODs). Such a distribution of metal-binding residues in sectors was also consistent with the taxonomic levels of the origins of these SODs, implying the functional separation of residues responsible for catalysis and structural stabilization. This study will help to explain that mutations apart from the active site can also exert pronounced effects on metal ion binding and provide more targets for engineering the biochemical properties of SODs.