Intragenomic polymorphisms among high-copy loci: a genus-wide study of nuclear ribosomal DNA in Asclepias (Apocynaceae)

Sampling and sequencing

One hundred twenty-five individuals representing 90 Asclepias species and subspecies were sampled (Table 1) and sequencing libraries were produced as described in Straub et al. (2012). Two individuals of putatively hybrid origin were included: A. albicans × subulata and A. speciosa × syriaca. These individuals were collected from wild populations and identified as hybrids through expression of intermediate morphological characteristics (M Fishbein, 1996, 1998, unpublished data; see also Fishbein et al., 2011). Samples were multiplexed in approximately equimolar ratios, with up to 21 individuals per lane, and sequenced with 80 bp single-end reads on an Illumina GAIIx instrument (Illumina, San Diego California, USA). Asclepias subverticillata was multiplexed in a lane with 32 samples and sequenced with 101 bp paired-end reads on an Illumina HiSeq 2000 instrument, with reads analyzed as though they were single-end. One individual of A. syriaca was sequenced at higher coverage: this individual was sequenced in a single lane on an Illumina GAIIx with 40 bp single-end reads (Straub et al., 2011). To allow more efficient assembly downstream, read pools were filtered to remove plastid reads (using the custom script sort_fastq_v1.pl modified to retain Ns; Knaus, 2010).

An A. syriaca haploid genome size estimate of 420 Mbp and a nrDNA copy number estimate of 960 were used for estimates of sequencing depth and for comparisons with other organisms. These estimates are modified from Straub et al. (2011), where an incorrect estimate of the average A. syriaca 2C value led to a haploid genome size estimate of 820 Mbp and a nrDNA copy number estimate of 1,845. The current values are based on 2C estimates from Bai et al. (2012) and Bainard et al. (2012).

Polymorphism quantification

The method for determining polymorphisms present among nrDNA copies within an individual while retaining information about position homology across a group of distantly related individuals included four general steps, detailed below: (1) A sequence was selected to serve as a reference for the whole group; (2) A consensus sequence was obtained for each individual taxon and aligned against the group reference, allowing tailored read-mapping for each individual while associating positions along the individual consensus with their homologous positions in the group reference; (3) Reads for each individual were mapped onto that individual’s consensus sequence; (4) At each position the reads differing from the individual consensus were tallied.

RNA structure determination

The secondary structure of each subunit and spacer region was predicted for the A. syriaca reference from the minimum free energy structure found by the program RNAfold for the 18S, ITS1, and ITS2 regions and RNAcofold for the 5.8S + 26S regions (Lorenz et al., 2011). Program default model parameters were used (37 °C, unpaired bases can participate in up to one dangling end). Positions along the cistron were then categorized as either paired (stems) or unpaired (loops). Predicted structures are provided in File S2.

Polymorphisms within the cistron

The effects of cistron position (subunit or spacer) and secondary structure position (stem or loop) on the likelihood of a position being polymorphic or highly polymorphic were assessed in two ways. Differences in the likelihood that at least one individual was polymorphic at a position (e.g., positions coded as either polymorphic or invariant) were assessed via a two factor multiple logistic regression, as implemented in R ver. 3.1.0 using the MASS ver. 7.3.33 package (Venables & Ripley, 2002; R Core Team, 2014). Differences in the abundance of polymorphic individuals at a position were assessed using square-root transformed data with a two-way ANOVA and type III sum of squares for unbalanced design, as implemented using the car ver. 2.0.20 package (Fox & Weisberg, 2011). The ANOVA analysis is not reported for the highly polymorphic individuals because the data diverge substantially from assumptions of a normal distribution.

Phylogenetic context

A maximum likelihood estimate of phylogenetic relationships within Asclepias was produced by Fishbein et al. (2011, Fig. 2, TreeBase #27576). This tree was pruned to match the sampling in this study, and counts of polymorphic positions were recorded for each taxon. Taxa sampled in this study, but not present in Fishbein et al. (2011), were omitted from further analyses. Counts were averaged for taxa with multiple individuals in this study, but sampled only once in Fishbein et al. (2011). Ancestral states for the number of polymorphic positions in the rDNA cistron were reconstructed using squared-change parsimony in the Mesquite phylogenetic suite ver. 2.75 + build 573 (Maddison & Maddison, 2011; Maddison, Maddison & Midford, 2011).

The phylogenetic signal in the distribution of the number polymorphic positions was tested in two ways. In the first method, the total length of the tree (parsimony steps) in terms of changes in number of polymorphic positions was compared to a distribution of tree lengths created by 10,000 permutations of polymorphic positions across tips. In the second method, a likelihood ratio test (LRT) was performed between models where character evolution followed a Brownian motion model across the tree. In the first model, the parameter lambda (describing how well the phylogeny correctly predicts the covariance among taxa for a trait) was found that maximized the model’s likelihood (Pagel, 1999). This was compared to the likelihood found when lambda was held at zero, representing phylogenetic independence among species for that trait. Parsimony permutations were performed in Mesquite (Maddison & Maddison, 2011; Maddison, Maddison & Midford, 2011), and likelihood ratio tests were performed in R with the phytools ver. 0.4.05 and ape ver. 3.1.2 packages (Paradis, Claude & Strimmer, 2004; Revell, 2012; R Core Team, 2014).

To determine if any clades held intragenomic polymorphism frequencies that were significantly high or low relative to the rest of the phylogeny, polymorphism rates were simulated along the tree, and the true polymorphism counts compared to the distribution of simulated counts (Garland et al., 1993). Two tips in the tree separated by zero branch length (A. asperula ssp. asperula and ssp. capricornu) were collapsed, and the average of the polymorphism counts for the four sampled individuals of the species used for the new tip. A model of trait evolution under Brownian motion, using the lambda parameter estimated from the LRT above, was found from 10,000 random starting points using the fitContinuous function from the geiger ver. 2.0.3 R package (Harmon et al., 2008). This model was used to simulate polymorphism counts across the phylogeny 10,000 times, with lower and upper bounds of 0 and infinity, respectively, using the fastBM function in phytools ver. 0.4.31 (Revell, 2012). For each node, the ancestral state was estimated for the true data and the simulated data using squared change parsimony as implemented in Mesquite (Maddison & Maddison, 2011; Maddison, Maddison & Midford, 2011).

Intragenomic polymorphism

All individuals were polymorphic at several positions homologous with the A. syriaca reference (Table 1). The number of polymorphic positions ranged from 23 (A. verticillata, 0.41% of sequenced bases) to 882 (A. solanoana, 15.13%), with a mean of 333 (5.77%, Fig. 1). A very high percentage (91%) of positions in the A. syriaca reference were polymorphic in at least one individual (Fig. 2A).

Positions were significantly more likely to be polymorphic if they were in a spacer region (ITS1, ITS2; Fig. 2A) or stem (Fig. 3A). This is true both when considering the number of individuals polymorphic at a position (Table 2), and whether any sample was polymorphic at that position (Table 3A).

The number of highly polymorphic positions ranged from 2 (A. arenaria, 0.03%) to 111 (A. boliviensis, 1.91%) with a mean of 28 (0.50%). Positions highly polymorphic in more than 10 individuals were found in the 18S, ITS1, ITS2, and 26S regions (Fig. 2B). The most polymorphic position was 4,172 (using the A. syriaca reference), in the 26S region, which was highly polymorphic in 29 individuals. Highly polymorphic positions were dramatically less frequent in the subunit regions (18S = 19%, 5.8S = 26%, 26S = 31%) than in the spacer regions (ITS1 = 61%, ITS2 = 57%; Table 3B, Fig. 2B). Positions in secondary structure stems were moderately more likely to be highly polymorphic in at least one individual than loop positions (Table 3B, Fig. 3B).

Relaxed read mapping

Allowing more mismatches when mapping reads to their consensus nrDNA sequence increased the polymorphic sites counted within individuals by an average of 15.5% (Fig. S2). Two samples had no change in their polymorphism abundance, and two had a decrease (i.e., newly mapped reads at a previously polymorphic position matched the consensus, thereby dropping the polymorphic reads below 2%). The increase in polymorphism abundance under relaxed read mapping may indicate that the standard read mapping parameters are too conservative and that some truly polymorphic sites are excluded. However, because standard read mapping is more likely to exclude reads containing sequencing errors, and because there is a strong linear correlation (R² = 0.97) between polymorphism abundance under the two mapping schemes, remaining analyses will only consider results from the standard read mapping.

Samples contained a mean of 7.6 and a median of 5 polymorphic indels when mapped reads were allowed to contain insertions or deletions of up to 5 bp relative to the individual consensus sequence. Eight individuals exhibited no polymorphic indels, including a sample of A. tuberosa ssp. rolfsii (Lynch 12526 [OKLA]). However, a sample of A. tuberosa ssp. interior (Fishbein 2816 [OKLA]) contained the most polymorphic indels at 51. Intragenomic indel abundance was positively correlated with SNP abundance (R² = 0.35, Fig. S3). Due to the generally low level of intragenomic indel polymorphisms (78% of samples contained 10 or fewer), remaining analyses only consider results from intragenomic polymorphic SNPs. Polymorphic indels and SNP counts under relaxed mapping for each sample are provided in File S3.

Phylogenetic signal

The number of polymorphic base pair positions exhibited strong phylogenetic signal across Asclepias under both the permutation test (P < 0.0012) and the likelihood ratio test (estimated lambda = 0.51, P = 0.0067; Table 4; Fig. 4). This signal remained even after the ITS regions were removed from the dataset (permutation test P < 0.0200, lambda = 0.45, LRT P = 0.0112). Highly polymorphic base pair abundance, however, was not significantly influenced by phylogenetic history (Table 4; Fig. 5), even when considering only the subunit regions.

Of the 53 resolved clades in the phylogeny, none showed polymorphism values more extreme than expected under a Brownian motion model after correcting for multiple comparisons. The most extreme clade was that formed by A. hypoleuca and A. otarioides, which had an ancestral number of polymorphic positions exceeded by 79 of the 10,000 simulations (when excluding the spacer regions this node was exceeded by 60 simulations). The following most extreme nodes were those ancestral to A. rosea and A. lemmonii with more polymorphic positions than all but 130 (79) simulations, and ancestral to A. boliviensis and A. mellodora with more polymorphic positions than all but 246 (150) simulations. Under a Bonferroni correction a clade would require 9 or fewer simulations more extreme than the observed value (α = 0.05) to reject a hypothesis of no divergence from Brownian motion.

Despite the strong phylogenetic signal of polymorphism abundance across Asclepias, counts among samples within species (for those species with multiple samples) exhibited variability. Some samples of the same species had very similar polymorphism counts (e.g., the two A. jaliscana individuals contained 140 and 165 intragenomic polymorphisms), while others differed dramatically (e.g., A. macrosperma individuals contained 76 and 391).

Identification of mixed ancestry

The number of polymorphic sites for hybrid individuals, 299 for A. albicans × subulata and 161 for A. speciosa × syriaca, are less than the mean number of 333 polymorphic sites; and the number of highly polymorphic sites, 12 and 37, are less than or greater than the mean of 28. Of those positions that are highly po lymorphic, 4 of 37 in A. speciosa × syriaca have a minor allele frequency of 0.3 or higher, while none of the positions in A. albicans × subulata have a minor allele frequency above 0.2.

Polymorphism patterns across the nrDNA cistron

Spacer regions (ITS1, ITS2) had higher frequencies of polymorphic positions than subunit regions (18S, 5.8S, 26S; Fig. 2). However, positions with low polymorphism frequencies are distributed much more evenly across the nrDNA cistron (Fig. 2A) than highly polymorphic positions, which show strong differentiation between the subunit and spacer regions (Fig. 2B). The lower frequency of highly polymorphic positions within the subunit regions suggests that these regions are under selection to remain homogenous within individual genomes. The lower difference in low polymorphism frequency between subunit and spacer regions suggests that this selection pressure is positively correlated with the proportion of nrDNA copies that differ from the majority. These findings contrast with those reported for nematodes (Bik et al., 2013), where the subunit regions had much higher levels of polymorphism abundance than the spacer regions.

Positions in stem regions were more likely to be polymorphic than loop positions (Fig. 3). This was strongly significant for all polymorphic positions (Tables 2 and 3A) and moderately significant for highly polymorphic positions (Table 3B). This would seem to contradict the hypothesis that stem sites in general should be more highly conserved in order to maintain a functional RNA secondary structure. Indeed, this finding agrees with those from Rzhetsky (1995), who not only found that trees estimated from stem regions contained longer branch lengths than those from loop regions, but that those stem sites least likely to affect secondary structure tended to be less conserved. Loop sites, on the other hand, contain a large proportion of the sites critical to ribosomal function (Rzhetsky, 1995), and may be under stronger stabilizing selection than stem sites.

Among the subunit regions, the 26S region harbored the highest frequency of highly polymorphic positions. This agrees with Stage & Eickbush (2007) who showed the Drosophila 28S region to have a higher mutation rate than the other subunit regions. However, in that study the 28S region also had a lower frequency of polymorphic base pairs. Stage & Eickbush (2007) hypothesize that this is due to the action of two retrotransposable elements found in many Drosophila 28S copies, with instances of aborted insertions causing the cell to repair the region using a nearby template and thereby homogenize the copies. They predict that levels of polymorphism will be higher in the 28S region of organisms lacking the retrotransposable elements, as seen here in the homologous angiosperm 26S region.

Within the subunit regions, highly polymorphic positions may be more common within expansion regions of the RNA gene (less conserved regions that tend to incorporate sequence insertions without affecting functionality; Clark et al., 1984). This is strongly implied by the recovery of clusters of highly polymorphic sites near A. syriaca positions 4,440 and 5,020, which are directly within the 26S expansion regions seven and eight, respectively (Kolosha & Fodor, 1990; Fig. 1). This agrees with results found in Drosophila (Stage & Eickbush, 2007). However, this may not be true for all highly polymorphic sites, as the high peak at position 4,172 is not within an expansion region.

Phylogenetic context

Mapping the number of polymorphic positions onto the phylogeny of Asclepias demonstrates strong and significant phylogenetic signal when counting all polymorphisms (Fig. 4), but this signal is not significant when only counting highly polymorphic positions (Fig. 5). While this study demonstrates phylogenetic signal in polymorphism abundance, it remains unknown whether abundance within a lineage is influenced by selection or purely neutral causes. Different demographic histories across clades could create this effect via neutral causes, while selection against organisms that retain too many variant nrDNA copies could be variable across the genus. Polymorphisms present in a high proportion of nrDNA copies may be uniformly selected against, while low polymorphism sites are tolerated at varying levels. This could explain the lack of phylogenetic signal at the highly polymorphic sites, and its presence when including all polymorphic sites.

These results have implications for phylogenetic inference based on nrDNA data. As previously cautioned, it cannot be assumed that all copies of nrDNA are identical within a genome. This is especially true for the spacer regions, but also for the subunit regions (Álvarez & Wendel, 2003). The discovery of phylogenetic signal in intragenomic polymorphism abundance demonstrates that those positions likely to lead to ambiguities are not distributed evenly across the phylogeny (Fig. 4). In addition to topology, variable polymorphism rates among lineages may affect the estimation of branch lengths on a tree. This may become especially problematic when nrDNA is used to date phylogenies. Recently developed methods of phylogenetic inference incorporating information about polymorphic positions may be able to alleviate difficulties in tree building caused by uneven levels of polymorphism abundance (Potts, Hedderson & Grimm, 2014).

Identification of mixed ancestry

Characterization of intragenomic nrDNA polymorphisms may allow for the identification of hybrid offspring between parents with differing nrDNA sequences (Zimmer, Jupe & Walbot, 1988), and may be able to provide an estimate of the number of generations since hybridization. This is especially true for early generation hybrids, as nrDNA homogenization can occur in a small number of generations (Kovarik et al., 2005), and mosaic Saccharomyces genomes have been identified based on intragenomic polymorphism abundance (West et al., 2014). Neither of the polymorphism profiles for the two wild-collected putative hybrid individuals in this study strongly indicate that they are early generation hybrids. This is in contrast to evidence from nuclear gene sequences that show heterozygosity consistent with inheritance of divergent alleles from the putative parents in the A. albicans × A. subulata hybrid (B Haack and M Fishbein, 2013, unpublished data). Detection of mixed ancestry may be hampered in this case by a lack of fixed differences between parental haplotypes. Consensus sequences for the A. albicans and A. subulata individuals in this study show no fixed differences, while one difference is found between the A. syriaca and A. speciosa consensus sequences. Unexpectedly, this position is monomorphic in the hybrid, while the positions with minor allele frequencies >0.3 are in positions not shown to differ between A. syriaca and A. speciosa. This likely indicates ancestry from an unsampled nrDNA haplotype.

Individuals of the same species with dramatically different polymorphism profiles may indicate the presence of cryptic diversity. Differing demographic histories between populations within a species may lead to individuals from those populations possessing different levels of intragenomic polymorphisms. As with the interspecific hybrids discussed above, early generation hybrids between populations within a species could also possess inflated levels of intragenomic polymorphisms.