We present a new method for detecting coevolving sites in molecules. The method relies on a set of aligned sequences(nucleic acid or protein) and uses Markov models of evolution to map the substitutions that occurred at each site onto thebranches of the underlying phylogenetic tree. This mapping takes into account the uncertainty over ancestral states andamong-site rate variation. We then build, for each site, a ‘‘substitution vector’’ containing the posterior estimates of thenumber of substitutions in each branch. The amount of coevolution for a pair of sites is then measured as the Pearsoncorrelation coefficient between the two corresponding substitution vectors and compared to the expectation under the nullhypothesis of independence. We applied the method to a 79-species bacterial ribosomal RNA data set, for which extensivestructural characterization has been done over the last 30 years. More than 95% of the intramolecular predicted pairs of sitescorrespond to known interacting site pairs.