High-Level Molecular Diversity of Copper-Zinc Superoxide Dismutase Genes among and within Species of Arbuscular Mycorrhizal Fungi

G. intraradices single-spore isolates were established from a population located in Switzerland. The 17 isolates examined originate from a small field site (90 m by 110 m) located at Hausweid, Tänikon. For details on spatial arrangement see Koch et al. (25). All isolates of G. intraradices (Tänikon population and within-population comparisons), isolate DAOM 181602 (interspecific comparisons), Glomus sp. (strain MUCL 43196), Glomus proliferum (MUCL 41827), Glomus geosporum, and Gigaspora rosea were clonally subcultured every 15 weeks in two-compartment plates with root-inducing transfer DNA-transformed carrot roots (29). Freshly isolated hyphae and spores of each isolate and species were separately dried overnight at 48°C and ground into a fine powder using a Retsch MM300 mixer mill (Retsch, GmbH). The DNA was extracted using a modified version of the Cenis method for fungal DNA extraction with an additional step of a 1:1 dilution with a solution of 24:1 of chloroform-isoamyl alcohol before the final precipitation to remove remaining impurities (9).

A total of four specific and slightly degenerate primers were designed based on the published Sod1 sequence of Gigaspora margarita (27). Primer sequences were the following: 5′-AAAGCTATTGCAGTTCTTAG-3′ and 5′-CCCAATAACACCACAGGCAA-3′; 5′-GCAGTTITTTASWCCTGATAA-3′ and 5′-TAACACCACAGGCAACACGA-3′.

DNA amplifications yielded one single fragment of the expected size with each pair of primers tested. The PCR products were separated on agarose gel and purified using a MinElute Gel Extraction Kit (Qiagen). The fragments were cloned with a TOPO TA Cloning Kit (Invitrogen, Inc.) following the manufacturers instructions. DNA sequences were determined using a version 3.1 Terminator cycler sequencing kit and separated on an ABI Prism 3100 genetic analyzer (Applied Biosystems).

Sequences were analyzed and aligned using the Vector NTI package (Informax, Oxford, United Kingdom). Chromatograms were carefully examined by eye. A phylogenetic framework was reconstructed based on all Sod1 genes and internal transcribed spacer (ITS) sequences using the neighbor-joining and minimum-evolution methods implemented in MEGA, version 3.1 (26). The maximum-likelihood tree was computed under the WAG+G model of evolution, chosen according to the Akaike information criterion (as computed with ProtTest, version 1.4 [1]), using Phyml, version 3.0 19). Bootstrap support was calculated from 1,000 replicates. Bayesian posterior probabilities were calculated using MrBayes, version 3.1.2, (22) with priors, chain number, and temperature set to default values; the aamodelpr parameter was fixed to particular protein models of evolution chosen according to the Akaike information criterion. Two parallel Markov chains were run for 2 × 106 generations; every 100th tree was sampled, and the first 5 × 105 generations were omitted from topology and probability reconstruction.

Codon-based likelihood ratio tests (LRTs), implemented in the PAML package, were used to detect potential signatures of selection in our data. Comparing the number of nonsynonymous to synonymous substitutions per site (dN/dS) is typically used to detect selection at the molecular level. Positive selection is inferred by a dN/dS ratio (ω) exceeding 1, while purifying selection is characterized by a ω of <1, and neutrality is indicated by a ω of 1. To test hypotheses of molecular evolution, we used models that allow different ω rates at different sites: the null model M1a (assuming two site classes: sites under purifying selection and neutrally evolving sites) and the alternative model M2a (adding a third site class: sites under positive selection) (30, 38). Comparison of the log-likelihood values under the two models reveals whether the model that allows for positively selected sites fits the data set better than the null model. The LRT compares twice the difference of the log-likelihood to a χ² distribution with a degree of freedom equal to the difference of free parameters in the corresponding models. Naïve and empirical Bayes calculations were used to analyze in more detail the site class of positive selection to identify particular sites under positive selection. The codon-based likelihood ratio tests have been performed among the six AMF species using the phylogenetic relationships represented in Fig. S1A in the supplemental material. The G. intraradices sequence used in the species data set corresponds to the isolate DAOM 181602. Rearranging the positions of G. proliferum, G. geosporum, and Glomus sp. within the AMF phylogeny did not affect the results reported in this study. The population data sets have been analyzed using the phylogenetic relationships shown in Fig. S1B in the supplemental material. The presence of recombination has been investigated using the program RDP2 (28) and default parameters. Recombination events have been investigated using an alignment of nucleotide sequences of the Sod1 gene from the totality of AMF species and isolates.

Real-time PCR was performed on three G. intraradices isolates (A2, A4, and C2) to compare relative copy numbers of the Sod1 genes. Using the Sod1 sequences we obtained from the 14 G. intraradices isolates, we designed primers that amplify a 100-bp fragment and that annealed to a region conserved among all variants. The primers and probe were as follows: forward primer (5′ to 3′), CGA TTA CAG AGT TGG GAC CAC TTA; reverse primer (5′ to 3′), CGT TGG CGA TTT GGG TAA TG; probe FAM (6-carboxyfluorescein), (5′ to 3′): AGC AAC TTT ACC ATC CGG GCC AGC. The probe was labeled with FAM at the 5′ end and with Black Hole quencher 1 at the 3′ end. In the real-time PCR, the fluorescence of the probe was measured at each cycle at the annealing phase of the reaction. Real-time PCR amplification with the FAM-labeled probe was performed in 40 μl containing 1× qPCR Mastermix (Eurogentec), a 0.5 μM concentration of each primer, 0.25 μM probe, and 10 μl of DNA at different concentrations. The thermal cycling conditions comprised an initial step of 50°C for 2 min and 94°C for 10 min, followed by 45 cycles of 94°C for 15 s and 55°C for 1 min. Fluorescence data were collected using an ABI Prism 7000 Sequence Detection System (SDS; Applied Biosystems). The SDS software then generated each real-time PCR profile after multicomponent analysis by plotting the log of the change in fluorescence (delta Rn, where Rn is the normalized reporter signal) versus cycle numbers. The cycle threshold (C_T) was determined by the SDS software as the fractional cycle number. This indicates where the fluorescence crosses an arbitrary threshold intersecting the signal curves in their exponential phases. In each experiment, twofold serial dilutions of G. intraradices genomic DNA from the different isolates (ranging between 5,000 and 39.06 pg) were included to generate C_T values. Two independent experiments were performed, each with eight concentrations of the genomic DNA. In half of the experiments, the DNA of the G. intraradices isolates was obtained from two independent extractions, thus allowing the possibility of checking whether a bias occurred in the C_Ts because of DNA extraction procedures. The C_T values shown in Fig. 3 represent averages between the two independent experiments. All replicates offered very similar results and produced standard curves with regression coefficients (R²) of >0.99. The relative copy number of the Sod1 gene in each of the G. intraradices isolates could then be calculated easily by comparing the C_T values of each of the isolates for a given amount of genomic DNA. The slopes resulting from the analyses of C_T values relative to the log of the amount of genomic DNA were consistent and reproducible among isolates, and the efficiency of the quantitative PCR from genomic DNA was found to be more than 95%.

AMF sequences obtained in this study have been deposited in GenBank under accession numbers EU232637 to EU232655.

To better understand the molecular evolution of the Sod1 genes in AMF, we sequenced 360 bp of the locus (almost the entire open reading frame) from five different AMF species (in two genera), as well as the genes from 14 different G. intraradices isolates harvested from the same field. All the sequences obtained showed a high degree of homology to the deposited sequences encoding fungal Cu-Zn superoxide dismutases according to a tBlastx search procedure (2). A phylogeny reconstructed using amino acid sequences from a variety of fungal Cu-Zn superoxide dismutase (Sod1) sequences, including the AMF sequences analyzed in this study, revealed the monophyly of all major fungal phyla, including the Glomeromycota (Fig. 1). The relationships we have identified between the fungal phyla are consistent with the most recent phylogenetic studies of the fungal kingdom (23). Interestingly, all the phylogenetic models used resulted in a basal positioning of the chytrids and the zygomycetes within the fungal kingdom; however, their monophyletic origin was not supported by maximum-likelihood and Bayesian methodologies. Overall, statistical support for the major fungal clades was relatively low, a likely consequence of high amino acid sequence diversity within and across the major fungal phyla and the small size of the studied locus.

Amino acid comparisons among all the AMF species and isolates revealed a substantially low average sequence similarity of 73% among the analyzed species. This suggested that these sequences might have evolved either under neutrality or positive selective constraints for some time. High sequence variation was also observed among 14 G. intraradices isolates harvested from one field. A total of five different sequences of Sod1 were obtained from 14 isolates of G. intraradices. Sequences from three of the isolates (A2, A4, and C2) accounted for most of the sequence variation, suggesting that the Sod1 locus diverged significantly in these isolates, a feature that is unexpected for individuals from a single population (Table 1).

This surprisingly high molecular diversity warranted independent controls to ensure that it was not caused by an erroneous species assignment and to possibly rule out paralogous origin of the Sod1 sequences from our population. These controls were performed as follows. All isolates of G. intraradices used in this study were identified through their ITS sequences, and their monophyly was confirmed by using the closest related species, Glomus diaphanum and Glomus fasciculatum (Fig. 2) (15). In parallel, the putative orthology of the Sod1 sequences and the absence of intraindividual variation were determined using two approaches. First, we screened and sequenced 32 clones of the Sod1 gene in every single isolate used in this study, using two independent sets of primers (see Materials and Methods); no sequence variation was found within isolates. Based on the sequencing results, we were very confident that the molecular diversity we identified was representative of variation among isolates and that only a single sequence variant was present in each isolate. The absence of intraindividual variation also suggested that the divergent Sod1 sequences did not arise through gene duplication (i.e., absence of paralogy). However, as cloning techniques are prone to randomly missing certain sequences and do not provide information about Sod1 copy numbers, we performed a real-time PCR experiment using an additional set of primers (see Materials and Methods) targeting a highly conserved region of the Sod1 gene in three of the most divergent isolates used in this study (Fig. 3). A significant difference in standard curves (one C_T difference corresponding to a twofold increase in gene copy number) would have suggested a difference in copy number for the Sod1 gene among the isolates (11, 13). The real-time PCR procedure showed very similar standard curves for all three isolates, suggesting that all G. intraradices isolates share the same copy number for the Sod1 gene.

Considering the amount of sequence variation we found among AMF species for the Sod1 locus, we tested whether this locus might have evolved either under neutrality or under the regime of positive selection. To test this, we analyzed the Sod1 sequences from the six different AMF species analyzed in this study using different models of codon evolution implemented in the PAML software package (30).

The results from the maximum-likelihood analysis of codon evolution showed that positive selection potentially acted at specific codons throughout the evolution of these AMF species. Indeed, all models that allowed for sites with a ω (dN/dS) greater than 1 (M1, M2, and M8), that is, the models of positive selection, fitted the Sod1 data significantly better than the corresponding neutral models (M0, M3, and M7) (Table 2). More specifically, the ω estimates, averaged over all sites and all branches, ranged from 0.5 to 0.61 among selection models. Models M1a and M7 gave lower estimates as they do not account for sites under positive selection. An LRT comparing model M1a (nearly neutral) to M2a (positive) indicated that the selection models fitted better (P = 0.009) (Table 3). The LRT comparing the model M7 (β, neutral) with M8 (β, selection) also rejected the null model (P = 0.004). Most sites in the fragment examined have apparently evolved either under purifying selection (M2a, 69%; M8, 92%) or neutrality (M2a, 26%; M8, 8%) (Table 2). The discrete model associated 10% of the sites under positive selection. Similarly, M2a and M8 indicated positive selection at a small subset of sites (4.7% and 7.7%) located within the fragment analyzed (Table 2). In more detailed analyses, a Bayesian approach identified seven sites affected by positive selection within the AMF lineage under model M8 and five under model M2a (Table 3). Overall, five amino acid sites were identified as potentially affected by positive selection when all sites that allowed for the presence of positive selection were used in the analysis (5G, 55F, 58T, 79P, and 88D). Sites 79P and 93D were predicted to be subject to positive selection at the 95% confidence level using all models.

In the Cu-Zn superoxide dismutase protein, the binding of copper and zinc ions requires the presence of histidine amino acids at specific conserved sites, and, as expected, these motifs did not show significant signs of positive selection. Apparently, the motifs potentially subject to positive selection are spread throughout the coding sequence, and none have been identified within the conserved internal catalytic cleft. Although the posterior probability for these mentioned sites to be under positive selective constraints was relatively high, we could not match any of these with a known function for the Sod1 protein. The significance of our tests might denote a functional importance for this rapid variation at the amino acid level, especially for the amino acids 55F and 58T, as these positions are located, respectively, 1 and 4 amino acids away from the to the catalytic cleft and have a 93% posterior probability (under model M3) of being shaped by events of positive selection.

Likelihood methods for detecting positive selection can be inaccurate if significant recombination has occurred among the sequences (4). The G. intraradices isolates analyzed in this study have been the subject of a multilocus genotyping study based on newly developed simple sequence repeat and nuclear gene intron markers (15). No events of recombination have been found between the isolates A2, A4, C2, and C5 using sequence data that was generated from these markers (14a). These isolates are the same ones that showed the highest molecular divergence in our study. To further test whether recombination played a significant role in shaping the diversity of the Sod1 sequences identified in our population, we performed an analysis to detect possible recombination events within the G. intraradices population (28). This analysis failed to identify any recombination events (data not shown), suggesting that other evolutionary forces played an important role in shaping the molecular diversity among the sequences.

To identify other possible evolutionary forces shaping the diversity at the G. intraradices Sod1 locus, we applied the same likelihood procedures reported in the previous section to the 14 G. intraradices isolates (which were initiated as single-spore cultures) harvested from one field. Consistent with what we identified among AMF species, we found that selection models applied to our population explained the data significantly better than neutral models (P < 0.05) (Tables 2 and 3). Although an LRT comparison between the discrete model M1a with M3 was only marginally significant (P = 0.09), the more stringent comparisons using M2, M8, and M8A (supplementary category ω fixed at 1) models were significant (P = 0.036) (Table 3). Importantly, models accounting for positively selected sites produced very similar parameter estimates and associated the same proportion (4.9%) of amino acids to the supplementary category of positively selected sites. Testing alternative tree topologies to take into account the phylogenetic incertitude did not alter the outcome of our analysis. Overall, sites at positions 6 and 18 (6G and 18G) were predicted to be under positive selection by all models at a high level of significance. Importantly, amino acids 6G and 58F were commonly identified as being under positive selection at both the population and species levels (Table 3).

AMF have sometimes been shown to harbor different alleles of one locus among the different genomes present in one isolate, but such intraindividual variability was not identified at the Sod1 locus. Although it could be important for the AMF to have copies with various functions for the Sod1 gene, the absence of such allelic variation is not surprising and has also been reported for other adaptively important genes, such as the plasma membrane-type II ATPases (13).

Rapid evolution of Sod1 has already been reported in pathogenic fungi (10). For example, the Sod1 enzymes from different varieties of the human pathogen Cryptococcus neoformans and Aspergillus sp. show substantial biochemical and physical differences (14, 21). In Aspergillus species, amino acid variation in the N-terminal domain of Cu-Zn dismutases has been shown to increase thermostability. This possibly affects their virulence as increased thermostability and activity have been reported in the most virulent strains with the potential to develop aspergillosis (21). While the role of the Sod1 genes in defense systems in plants, immune system in animals, and the pathogenesis of some fungi has been well characterized, their molecular variability has rarely been discussed in the context of adaptive evolution.

Our results are consistent with a rapid evolution of the Sod1 gene, which has been shown to be specifically upregulated during mycorrhizal symbiosis and upon exposure to root exudates. However, its specific role during the establishment of the mycorrhizal symbiosis is still unclear. For instance, it is not known whether its expression is directly related to the plant-fungus dialogue or whether it is in response to a change in environmental conditions dictated by the plant defense system (i.e., expression required to counteract the ROS released by the host). Rapid evolution of genes in a plant-fungal interaction might also be an indication of arms races, but such arms races are usually seen in host-pathogen interactions (6) and not expected to evolve in a mutualisitic symbiosis. Thus, an arms race scenario, although possible in the case of the Sod1 gene, seems unlikely considering the generalist nature of the AMF. From complementation studies in yeast, it has been demonstrated that the AMF Sod1 encodes an active polypeptide that effectively scavenges ROS (27). Moreover, this almost certainly occurs in AMF-colonized roots because suppression of Sod1 activity reduces build-up of peroxidase in the roots (27). It seems evident, therefore, that one important role for AMF Sod1 is the scavenging of ROS. This adaptive scenario implies that AMF are facing environmental situations similar to pathogens and, therefore, need to counteract the same plant defense barriers upon establishment of the mycorrhizal symbiosis. In this case, positive selection on the copper-zinc superoxide dismutase protein is better explained as a direct consequence of adaptation directed toward scavenging ROS produced by the host.

In parallel, Sod1 proteins have also been suggested as playing an important role in the fungal-plant molecular dialogue (18, 27) although conclusive evidence has not been provided for this. In this communication function of the Sod1 proteins, the rapid evolution of their amino acids could be related to the fact that tightly coevolved mutualisms are constantly exposed to nonmutualistic cheaters or parasites. As Sod1 plays a key role for pathogens to avoid plant immune systems (production of ROS), the AMF might need to slightly change their biochemical properties accordingly to avoid their proteins' being recognized as of pathogen origin. An alternate explanation could be related to the poor specificity of AMF species toward their host. High rates of protein evolution facilitating host interactions could render AMF species more broad in their host range or allow them to rapidly colonize hosts in the presence of competitors.

The molecular evidence we have gathered in this study points toward an orthologous origin of the Sod1 sequences we identified within the G. intraradices population. Although it is true that we cannot conclusively reject a paralogous origin for a few of these, the absence of AMF isolates and species harboring more than one Sod1 sequence and the impossibility to absolutely reject the monophyly of our G. intraradices sequences are both additional lines of evidence that make this latter evolutionary scenario far less parsimonious than the one suggesting an accumulation of mutations among orthologous sequences.

Intraspecific sequence polymorphism for a given gene has most often evolved as a consequence of random mutations accumulated by genetic drift. In this case, amino acid replacements are considered as neutral, with no effect or weak effects on protein function. An important departure from neutrality is observed when recombination produces molecular diversity or when positive selection fixes advantageous mutations. In the latter case, selective sweeps of favorable mutations lead to an excess of nonsynonymous over synonymous substitutions between species.

In this study, the molecular diversity identified among spores of G. intraradices is unlikely a consequence of recombination as other studies and our own analyses failed to identify recombination events. However, even if recombination did account for part of the intraspecific variation, this would not have significant consequences on our codon-based likelihood analyses for two independent reasons. First, rare events of recombination have been demonstrated as not interfering with the detection of positive selection (4). Second, LRT comparisons using the more realistic model M7 against M8 were shown to be robust to false positives due to recombination. This is particularly the case if the model M7 allows the positive selection pressure to vary between 0 and 1 (not accounting for positive selection) and if the alternative model M8 allows an additional discrete class with a dN/dS that could be estimated to be >1 (and thus account for positive selection) (4). In this study, the LRT comparison between these two models was the most significant one. This suggests that even under rare recombination in the population, diversifying selection is certainly a likely explanation for the diversity observed at the Sod1 locus.

Within a population, directional selection is thought to shift allele frequency in favor of a newly derived mutation through a fitness advantage to certain individuals over others. Signals of positive selection are, therefore, rarely detectable at the population level. However, a changing environment, with variable selective pressures over time (i.e., interactions with multiple symbiotic partners), could lead to balanced polymorphism that is maintained within a single population. This feature is known for genes involved in plant defense systems, mating systems in fungi and plants, or immune response in animals, where a high level of polymorphism is favored and maintained through selection (31). Thus, rapid protein evolution and balancing selection at the locus Sod1 may offer evolutionary advantages for some isolates compared to others in a heterogeneous environment, resulting in the maintenance of high protein divergence within the population.