Muons penetrate much further than X-rays, they do essentially zero damage, and they are provided for free by the cosmos.
Image credit: Science Source/Digital Globe.
Edited by David M. Hillis, The University of Texas at Austin, Austin, TX, and approved September 26, 2017 (received for review April 24, 2017)
Significance
Convergent and parallel evolution (homoplasy) is widespread in the tree of life and can obscure evidence about phylogenetic relationships. Homoplasy can be elevated in genomes because individual loci may have independent evolutionary histories different from the species history. We sequenced the genomes of five cacti, including the iconic saguaro of the Sonoran Desert and three other columnar cacti, to investigate whether previously uncharacterized features of genome evolution might explain long-standing challenges to understanding cactus phylogeny. We found that 60% of the amino acid sites in proteins exhibiting homoplasy do so because of conflicts between gene genealogies and species histories. This phenomenon, termed hemiplasy, is likely a consequence of the unusually long generation time of these cacti.
Abstract
Few clades of plants have proven as difficult to classify as cacti. One explanation may be an unusually high level of convergent and parallel evolution (homoplasy). To evaluate support for this phylogenetic hypothesis at the molecular level, we sequenced the genomes of four cacti in the especially problematic tribe Pachycereeae, which contains most of the large columnar cacti of Mexico and adjacent areas, including the iconic saguaro cactus (Carnegiea gigantea) of the Sonoran Desert. We assembled a high-coverage draft genome for saguaro and lower coverage genomes for three other genera of tribe Pachycereeae (Pachycereus, Lophocereus, and Stenocereus) and a more distant outgroup cactus, Pereskia. We used these to construct 4,436 orthologous gene alignments. Species tree inference consistently returned the same phylogeny, but gene tree discordance was high: 37% of gene trees having at least 90% bootstrap support conflicted with the species tree. Evidently, discordance is a product of long generation times and moderately large effective population sizes, leading to extensive incomplete lineage sorting (ILS). In the best supported gene trees, 58% of apparent homoplasy at amino sites in the species tree is due to gene tree-species tree discordance rather than parallel substitutions in the gene trees themselves, a phenomenon termed “hemiplasy.” The high rate of genomic hemiplasy may contribute to apparent parallelisms in phenotypic traits, which could confound understanding of species relationships and character evolution in cacti.
Cactaceae have undergone adaptive radiation on a continental scale in the Americas. Occurring from arid deserts to alpine steppes and tropical forests, they exhibit a remarkable diversity of growth forms, ranging from tiny, nearly subterranean “buttons” to giant columnar or candelabra forms, epiphytes, and leafy shrubs (1). Classification of the family’s 1,438 species (2) has been unusually fraught, with taxonomic treatments recognizing anywhere from 20 to 233 genera (1⇓⇓–4), and classification above the genus level being equally problematic (1, 5). In part this has been attributed to homoplasy (ref. 1, p. 18 and ref. 6, p. 45), the independent evolutionary origin of the same trait (7). For example, early taxonomists combined the giant columnar cacti of North and South America into one large genus Cereus Mill (3). Later split into numerous smaller genera, this association was sometimes maintained at a suprageneric rank (8), but molecular phylogenies clearly show North and South American Cereus are separate clades (6, 9).
Frequent convergence in growth habit may reflect design limitation (10) of the relatively simple cactus body plan of stem succulence, simplified branching, and loss or reduction of leaves (9) (ref. 11, p. 536). In cacti, convergent simplification via paedomorphosis (12), along with parallelisms among close relatives, has also obscured phylogeny (13, 14). This may help explain the finding that <50% of (nonmonotypic) cactus genera are monophyletic (15). Taxonomic impediments take on special significance in cacti because of the unusually high fraction of species of conservation concern (31% estimated in ref. 16). Despite its potential significance, homoplasy has been quantified in cacti for only a small number of phenotypic traits (14, 17). Here, we focus on the genomes of cacti and assess the degree to which apparent homoplasy among these species is elevated due to discordance between gene trees and the species tree.
An important but recently recognized contributor to molecular homoplasy is “hemiplasy” (18, 19), originally defined as the apparent multiple origin of a character state on the inferred species tree arising when an inferred gene tree (with no homoplasy) is discordant with that species tree (18⇓–⇓21) (Fig. 1). A slight generalization of this definition is needed to account for characters with homoplasy on both trees (Fig. 1, Lower Right): A character (site) exhibits hemiplasy if the inferred number of character state changes on the species tree is strictly greater than on the gene tree, which occurs only if the two trees are discordant (Materials and Methods).
Fig. 1.
Hemiplasy in amino acid alignments of cacti (genes annotated as Cgig1_gene##:site ##). Hemiplasy occurs when the inferred number of state changes of a trait is greater on the species tree than on the gene tree, which happens only if the two trees are discordant. Gene trees in blue are imbedded in black species trees. Yellow dot is the position of a state change inferred for an amino acid site for this gene tree. Black rectangles are locations where state changes would be inferred on the species tree (alternative equally optimal reconstructions are possible but the number of state changes is the same). Amino acids are indicated at leaf nodes. Three trees have the same number of state changes on species and gene trees, but at Lower Left a homology (one change) on the gene tree is seen as homoplasy (two changes) on the species tree because of discordance (i.e., hemiplasy). Cgig, C. gigantea; Lsch, L. schottii; Phum, P humboldtii; Ppri, P. pringlei; Sthu, S. thurberi.
Genome scale datasets revealed high levels of gene tree discordance (20, 22⇓–24), which suggests a potentially important role for hemiplasy as a contributor to overall molecular homoplasy. Discordance can arise because of incomplete lineage sorting (ILS) and gene flow, which depend, in turn, on demographic factors at the population level, phylogenetic history, divergence times, mutation rate, generation time, and the timing and taxa involved in introgression or hybridization (25). It can also arise from gene duplication and loss (26).
The cactus tribe Pachycereeae is a model for taxonomic and phylogenetic complexities in cacti (27, 28). Its ∼70 species, including most of the columnar cacti of Mexico and adjacent regions, have been dispersed among 6–23 genera in various taxonomic treatments (1, 5, 27). Despite early molecular support for its monophyly (11), broadly sampled recent molecular studies have led to abandonment of the tribe (9, 13, 15), or to recasting it with the informal name, “core Pachycereeae” (6). Notably, homoplasy in vegetative and floral traits has been cited as a factor contributing to its difficult taxonomic history (ref. 28, p. 1086 and ref. 29, p. 556).
To sample widely in Pachycereeae, we sequenced the genomes of four representatives of the two subtribes of Pachycereeae recognized in Gibson and Horak (30) (Table S1) and some earlier treatments. This included the iconic saguaro cactus of the Sonoran Desert (Carnegiea gigantea) to serve as a high-coverage reference assembly, and Pachycereus pringlei (cardón, sahueso), Lophocereus schottii (senita), and Stenocereus thurberi (organ pipe, pitaya), also of the Sonoran Desert, to lower coverage. Pereskia humboldtii, a leafy, Andean cactus, was included as an outgroup (31).
View this table:
Table S1.
Sample, voucher, and sequence library information for taxa in this paper
Results and Discussion
The assembly of saguaro’s 1.4-Gb genome (32) from short read libraries (Table S1) spanned 980 Mb, with a scaffold N50 of 61.5 kb (Table S2). Transcriptomes were also assembled (Table S3) and used with other evidence to annotate the saguaro genome. The saguaro genome contains 28,292 protein coding genes; 58% of the assembly consists of transposable elements and retroviral and repeated sequences (Tables S2 and S4). Assemblies from single libraries of the other four cacti were more fragmented (Table S5).
View this table:
Table S2.
C. gigantea (SGP5) assembly metrics (percent of assembled genome included where appropriate)
View this table:
Table S3.
Saguaro transcriptome sequence and assembly results
View this table:
Table S4.
Repeat annotation of the saguaro genome
View this table:
Table S5.
Comparison of assembly results for five cacti in this study
The high completeness of the gene space in the saguaro assembly (Table S2) let us construct 4,436 alignments of orthologous genes across the five cacti, each of which contained two or more exons for all five taxa, and contiguous introns. Based on gene tree confidence levels estimated from alignments containing all nucleotide positions, sets of alignments differing in phylogenetic robustness (“gene confidence sets”) were compiled for downstream analyses. For example, the 90% gene confidence set comprised 458 genes, with gene trees having maximum likelihood bootstrap support value above 90% for all clades (Table 1). These gene sets provide a control for the effect of weakly supported gene trees on gene tree discordance estimates (33).
View this table:
Table 1.
Gene tree discordance relative to the Pachycereeae species tree (Fig. 2B)
The most robust 90% gene confidence set implies levels of gene tree discordance between Pachycereeae genera (Fig. 2A) comparable to those seen in closely related short-lived angiosperm genera such as Solanum sect. Lycopersicon [time scale of 2 million years (MYR); ref. 24] and the Oryza AA genome clade (2.5 MYR: ref. 22). Species trees constructed using both gene tree-based (MP-EST; ref. 34) and alignment-based methods (BPP; ref. 35) were the same for all gene tree confidence sets and for three different partitioning schemes for the alignments by codon position and intron (Materials and Methods and Fig. 2B). An exhaustive search of tree space indicated a single optimal species tree under the MP-EST pseudolikelihood criterion, and three replicate runs of each MCMC chain in BPP from random starting trees generated this same species tree for all five gene confidence sets, indicating convergence. This optimal tree agreed with previous molecular phylogenetic analyses of mostly plastid genes (6, 28, 29). In the 80% and 90% gene confidence sets, 61% and 63% of gene trees, respectively, agreed with the species tree (Table 1). Estimates of species tree edge lengths for these gene confidence sets in BPP ranged from 1.11 to 1.24 for edge A and from 1.22 to 1.34 for edge B, in “coalescent time units” (= T/2Ne, where T is measured in generations), depending on data partition, with slightly lower values estimated by MP-EST (Table 1).
Fig. 2.
(A) Topological discordance in the 90% gene confidence set of 458 gene trees visualized in DensiTree (78). Blue trees are the most frequent, followed by red and shades of green. (B) Species tree inferred with BPP and MP-EST, which was the same for all gene sets, with divergence times estimated using the Neutral partition. Blue dashed line indicates reticulation, having inferred inheritance probability, γ, on the optimal PhyloNet analysis of the same set of gene trees. Green dashed line indicates position of reticulation on the next best inferred network (Fig. S1). The origin of this reticulation edge earlier than its endpoint implies the existence of an extinct or unsampled taxon.
Introgression may also contribute to gene tree discordance. Network reconstruction allowing both ILS and gene flow in PhyloNet (36) indicated substantial support for a network model with one reticulation vs. a tree model (ΔAIC = 17.0 and 28.6 for the 90% and 80% gene sets, respectively; ref. 37), implicating introgression into P. pringlei from either S. thurberi or some more closely related but unsampled or extinct taxon (Fig. 2B and Fig. S1). The best two networks were the same in both gene confidence sets (although the ranking reversed; Fig. S1). Moreover, the “major tree” within these networks (Materials and Methods) was the same as that inferred by species tree methods, and the inheritance probability estimate from the 90% gene confidence set, which indicates the level of introgression, was low (6.5% and 14.2% for the optimal and next best networks in the 90% gene confidence set; although higher in one of the two networks from the 80% gene set). The two optimal networks were consistently recovered in multiple searches from different starting networks, and in each gene set, the two results were substantially better than other suboptimal networks (Fig. S1).
Fig. S1.
PhyloNet (36) network reconstructions for the 90% (A) and 80% (B) gene confidence sets. Likelihood scores decrease from left to right. AIC scores are relative to the optimal networks at Left. Branch lengths are not scaled to sequence divergence, but the lesser of the two inheritance probabilities for reticulation edges is given next to edge in the format ##.#%. Parametric bootstrap probabilities are given next to edges for top scoring networks in A and B, as percentages in the format ##. The dotted box encloses the two best networks in each gene confidence set, and these are the same in A and B, but in reverse ranking. The remaining networks have substantially worse AIC scores relative to the optimal networks in each row.
Intergeneric hybrids are well known in cacti (1, 5), including rare hybrids between P. pringlei and Bergerocactus emoryi in Baja California in a narrow zone where the two species are sympatric (1, 38). Bergerocactus (not sampled here) is more closely related to Carnegiea, Lophocereus, and Pachycereus than to Stenocereus in the most complete cactus phylogeny to date (6). P. pringlei is possibly a recent autotetraploid (39), so postzygotic inviability and sterility barriers (40) would have had to have been overcome if the introgression took place following tetraploidy.
Owing to the limited level of introgression, we inferred demographic history using BPP assuming that ILS is the primary cause of discordance. The neutral mutation rate was estimated from the substitution rate in the “neutral” partition of the alignments (Table S6) assuming a root age of the tree at 26.88 MYR (Materials and Methods). Estimates ranged from μG = 5.98 to 6.07 × 10−8 per site for each generation across the five gene confidence sets. An alternative estimate based on pairwise synonymous Ks distances in CDS regions between Carnegiea and Pereskia from all 4,426 alignments was somewhat higher at μG = 8.75 × 10−8 per site per generation (Table S7). Per year rates are comparable to several angiosperm tree species (41⇓–43), but generation times in the cacti, ranging from 20 to 75 y (44⇓⇓–47), are 2–12 times longer. The BPP estimates of scaled mutation rate (θ = 4 NeμG), imply an Ne of 24,000–39,000 and 31,000–36,000 for edges A and B, respectively, over the 15 gene sets and partitions (Table S6). Combining the BPP estimates of θ with the pairwise Ks estimate of mutation rate produces ancestral Ne of ∼2/3 the BPP estimates. These ancestral Ne values are higher than estimates of recent Ne for some perennial angiosperms [Populus trichocarpa: ∼4,000−6,000 (48); Eucalyptus grandis: ∼11,000 (41); Amborella trichopoda: 5,000 (49)] but comparable to Populus balsamifera (44,000−59,000; ref. 50); and much less than Pinus taeda (560,000; ref. 51).
View this table:
Table S6.
Demographic inference based on BPP analysis
View this table:
Table S7.
Neutral mutation (substitution) rate analyses across cacti
The probability that a rooted gene tree with three taxa disagrees with the species tree that contains it is 2/3 e−T/2Ne, where T is the time in generations along the internal edge of the tree (52). This can be high if T is small, either because divergence time in years is small or generation time is large. Thus, the level of gene tree discordance we inferred in columnar cacti having Ne of 25–40,000, along edges with a duration of 2–3 million years would be far too high were it not for the long generation times of these plants, similar to findings in long-lived conifers (53).
Gene tree discordance has a strong impact on the distribution of protein sequence variation among these genera of cacti in generating hemiplasy. Taking the 80% and 90% gene confidence sets as samples of potentially phenotypically relevant amino acid sequence variation across these genomes, we partitioned the homoplastic amino acid sites in genes on the species tree into three parts: the fraction arising from homoplasy on concordant gene trees (Fig. 1, Upper Right), the fraction arising solely from gene trees with less homoplasy than the species tree but that are discordant with the species tree (hemiplasy: Fig. 1, Lower Left), and the small fraction that arises from homoplasy on discordant gene trees that is equally homoplastic on the species tree (which as yet has no term defined: Fig. 1, Lower Right). Hemiplasy accounts for 58–63% of all apparent homoplasy (Table 2).
View this table:
Table 2.
Homoplasy and hemiplasy in amino acid alignments of genes
As examples, consider two genes from the 90% gene confidence set having strongly discordant gene trees. Both play potential roles in the physiological adaptation of cacti to arid environments. The saguaro nuclear gene annotated as a chloroplastic NADP-dependent malate dehydrogenase (MDH: Cgig1_18427) is the only 1 of 11 MDH isoforms in the saguaro annotation that is NADP-dependent. It catalyzes the reduction of oxaloacetate to malate in the chloroplast (54) and appears to participate in the fixation of carbon dioxide under both light and dark conditions (55). Together with the cytosolic NAD-dependent MDH, these two isoforms may be primarily responsible for malate formation during dark CO2 fixation in crassulacean acid metabolism (CAM) plants (56). The MDH gene tree has a single replacement of an ancestral serine with an alanine, but the gene tree is discordant with the species tree, making the alanines appear to have evolved twice on the species tree, when they are in fact hemiplastic.
The gene annotated as a DNAJ JJJ1 homolog (Cgig1_00352) contains a DNAJ domain involved in heat shock protein interactions. Columnar cacti must be able to tolerate internal tissue temperatures that can exceed 50 °C in the summer (5). The gene tree is discordant with the species tree, and 19 of its 20 potentially informative amino acid replacements exhibit homoplasy on the species tree, but only one does on the gene tree, meaning 18 sites are hemiplastic. A search of the National Center for Biotechnology Information (NCBI) CDD database (57) shows that two of these hemiplastic sites are in the conserved DNAJ domain proper, although neither involves known HSP70 interaction sites.
Hemiplasy is not restricted to the coding sequence in this gene. The phylogeny of the 1,000 bp upstream of the start codon, which often contains proximal conserved regulatory elements in plants (58), is the same discordant tree as is found for the coding region. Of 22 potentially informative nucleotide sites there, 21 are homoplastic on the species tree, whereas only two are on the gene tree. Most homoplasy in this noncoding region on the species tree is thus also due to hemiplasy. All of these cases suggest caution in viewing multiple origins of the same amino acid or nucleotide in a species tree as possible evidence of convergent evolution (18).
These conclusions depend on the robustness of various model assumptions. The extent of gene tree discordance and hemiplasy was estimated from gene tree topologies inferred with an HKY substitution model in PAUP*. Estimates of tree topology are generally more robust than estimates of rates, especially at low sequence divergences (59⇓–61), where simple models sometimes even outperform more general ones (62). Mean pairwise nonsynonymous and synonymous distances between the most divergent Pachycereeae were only 1.3% and 4.1% respectively, and about one-half of the 4,436 sequence alignments had bootstrap values below 50%, a consequence of low sequence divergence and short length. Substitution model therefore likely had little impact on estimates of the prevalence of discordance and hemiplasy.
However, our explanation for these levels of hemiplasy in terms of generation time rests on estimates of divergence times and Ne obtained from BPP, which adds the assumptions of panmixia within and no gene flow between species, constant population sizes within tree edges, and no linkage (61). To examine the impact of BPP’s Jukes–Cantor clock model on its estimates of divergence times and effective population sizes, we evaluated more general models with subsets of the data using BEAST (63), finding that the parameter effect sizes were small (SI Materials and Methods and Table S8). Other assumptions were more difficult to test. Panmixia will not be strictly true even for outcrossing columnar cacti, but gene flow is likely high within three of the four Pachycereeae that are bat pollinated (e.g., S. thurberi; ref. 64). We did find limited gene flow between species in PhyloNet analyses, and this can lead to overestimates of population sizes and underestimates of divergence times when using the multispecies coalescent approach (65), but because the rates of ILS correlate with the product of population size and time, these biases potentially counteract each other. However, it is unlikely that population sizes have remained constant within species tree edges. Paleoclimatic and packrat midden evidence document the ebb and flow of Sonoran Desert plant communities during Pleistocene glacial and interglacial periods (66). BPP’s divergence time estimates are biased toward the recent in “more extreme bottleneck scenarios” of population history (67), but shorter edge durations in Pachycereeae would tend to make the observed level of ILS even higher than estimated. Finally, linkage is unlikely to be problematic. Our 80% gene confidence set has an average of 1.8 Mb between genes. Linkage disequilibrium in long-lived, outcrossing angiosperms typically decays in distances from a few kb to 50 kb (41, 48).
View this table:
Table S8.
Effect of model choice on parameter estimates in nine samples of 50 genes from 90% gene confidence set, assayed in BEAST 2 (63)
The phylogenomic history of saguaro and its relatives exhibits extensive gene tree discordance due to ILS and low rates of introgression. Comparable findings have been seen in rapid and/or very recent radiations (20, 22, 24), but in Pachycereeae, ILS acts at long time scales, between taxonomically divergent genera with very long generation times. A consequence of this discordance is elevated levels of apparent homoplasy in the species tree. The connection between apparent genomic homoplasy arising from ILS and apparent phenotypic homoplasy is probably strongest for traits with a simple genetic architecture (18), such as those involving the function or regulation of a single enzyme (21, 68). In plants, enzymes in biosynthetic pathways for floral pigments or defensive compounds are good candidates (69, 70). Notably, the taxonomic distribution of alkaloids, triterpenes, and sterols have played a role in the systematics of Pachycereeae (30). When hemiplasy arises because of introgression, genetic architecture may be less of an issue because multiple loci may undergo gene flow almost simultaneously (18). In plants, even complex traits with multiple components and fitness effects have been adaptively introgressed (71). Collectively, this body of evidence lends plausibility to the hypothesis that the phenotypic effects of genomic hemiplasy may have exacerbated the long-standing problem of inferring relationships in these charismatic cacti.
SI Materials and Methods
Sampling and Genome Sequencing.
Fresh material was obtained from an individual C. gigantea (Engelm.) Britton & Rose (“SGP5”) growing at the Tumamoc Hill Reserve in Tucson, AZ, and from a second plant purchased at a Tucson nursery (“SGP3”). Fresh material of L. schottii (Engelm.) Britton & Rose [=Pachycereus schottii (Engelm.) D. R. Hunt], P. pringlei (S. Watson) Britton & Rose, P. humboldtii Britton & Rose (=P. horrida DC), and S. thurberi (Engelm.) Buxb. were obtained from cultivated individual plants at the Desert Botanical Garden (Table S1). Genomic DNAs were extracted using a modified cetyltrimethylammonium bromide (CTAB) method (64).
For C. gigantea, paired end (PE) and mate pair (MP) libraries were constructed as described previously (82) (Table S1). The C. gigantea genome was sequenced in four different libraries (insert size ranging from 180 bp to 2.5 kb) on Illumina HiSeq (2 × 100-nt reads) and MiSeq (2 × 300-nt reads) sequencers. Single PE libraries were constructed for Lophocereus, Pachycereus, Pereskia, and Stenocereus and sequenced on Illumina NextSeq 500 (2 × 150 nt reads; Table S1). All reads were deposited in GenBank (Table S1).
Genome Assemblies.
Raw reads were trimmed using Sickle (83) for bases below QV20 and shorter than 20 bp. Genome size was estimated from k-mer frequency using the script KmerFreq_AR from the SOAPdenovo v2.01 suite (84). The C. gigantea data were error-corrected with BLESS v. 0.14 (85), and residual adapter bases were removed with PRINSEQ (86).
For all species, preprocessed reads were assembled de novo in scaffolds with Meraculous v. 2.0.4 (87) with the diploid option enabled. Out of the multiple iterations tested at different k-mer values, the most contiguous saguaro assembly was obtained with a k-mer value of 69; the genomes of P. pringlei, L. schottii, and S. thurberi were assembled with k-mer values of 55, 37, and 57, respectively. Gaps were closed with Platanus (88). Scaffolding was performed with BESST v. 1.2.3 (89) using the paired end and (when available) mate pair libraries; each scaffolding step was followed by a step of gap closing. Low confidence contigs or gaps were broken with REAPR v. 1.0.17 (90), producing the final assembly. Finally, we removed scaffolds and contigs shorter than 1 kb (saguaro) or 200 bp (all other species) that matched the ϕX174 genome, or that aligned to the saguaro chloroplast (82) for 90% of their length and had at least 90% similarity.
Blastx alignments against GenBank databases revealed significant hits only to other plant sequences. Completeness of the gene-space component of the genome assembly was assessed by running CEGMA (91) and BUSCO (92), using for the latter the Plantae database and tomato as the species for the Augustus gene prediction.
Annotation.
A saguaro repeat library was developed using RepeatExplorer (93) parsing the output as described previously (94). The annotation of the repeats in the saguaro assembly was obtained merging the output of RepeatMasker (www.repeatmasker.org/, v. 3.3.0) and Blaster (a component of the REPET package; ref. 95). Reconciliation of the two masked repeat sets was carried out using custom Perl scripts and formatted as gff3 files.
To identify noncoding RNAs, Infernal (96) was run using the Rfam library Rfam.cm.12.2. Hits with an e value higher than the threshold of 1e-5 were removed, as well as results with score below the family-specific gathering threshold. When overlapping loci were predicted, only the hit with the highest score was kept. Transfer RNAs were predicted using tRNAscan-SE (97) with default parameters.
We also obtained transcriptome data for saguaro using an RNA-seq approach. Samples were collected from SGP3 and from seedlings of the SGP5 specimen (Table S3). Libraries were constructed with the Ovation RNA-Seq System V2 (NuGEN) kit and sequenced on an Illumina instrument (Table S3). All reads were deposited in GenBank (Table S3).
The saguaro genome assembly was annotated using the MAKER pipeline (98). Several transcript and protein sequence datasets were identified to aid in gene prediction. Transcript unigenes were generated for the peyote cactus Lophophora williamsii from publicly available RNA-seq datasets (SRR1575212, SRR1575317, and SRR1575216). Saguaro RNA-seq reads were assembled de novo with Trinity r2013-02-25 (99), with the normalization option (Table S3). Assembled transcripts were merged in a unique EST file and clustered at 95% with Usearch v4.2.66 (100).
Protein sequences used for evidence during gene predictions included predicted proteins from Arabidopsis thaliana (TAIR10, ref. 101), Uniprot/Swissprot proteins (102), and GenBank protein sequences for the Caryophyllalean species Beta vulgaris and Spinacia oleracea. After repeat masking and transcript and protein alignment, high-quality saguaro transcript alignments were used to create an initial set of gene models that were used to train the Hidden Markov Model (HMM) for the SNAP gene prediction program (103). MAKER was run a second time, and genes were predicted using SNAP. High-quality SNAP gene predictions (annotation edit distances ≤ 0.2; ref. 104) were used to train the HMM for SNAP a second time, and MAKER was run again allowing the second SNAP HMM to predict genes. High-quality genes from the second SNAP run were then used to train the HMM for AUGUSTUS (105). MAKER was run a final time allowing genes to be predicted by both SNAP and AUGUSTUS. Alignments of single exon transcripts spanning less than 500 bp were excluded from the final MAKER run, and the “always_complete” and “keep_preds” options were set to 1. Predicted genes with 50 amino acids or less were discarded. Overall, 87% of the genes (25,862 models) had annotation edit distance between 0 and 0.5, with more than half (15,635) between 0 and 0.2.
Predicted proteins were analyzed with HMMER hmmscan to identify matching Pfam domains (106, 107). Gene models with support from transcript evidence, protein evidence, and/or matching Pfam domains were retained as an initial high-quality gene set that was further probed to identify and remove transposable element-related genes. Transcripts with homology to known transposable elements (TEs) and containing domains matching TE-related Pfam domains were excluded. The final set of 28,292 predicted genes were functionally annotated using the Trinotate pipeline (trinotate.github.io). Genes were deemed as transcribed if they had a Salmon v.0.6.1 (108) transcripts per million score equal or greater to 5, or if the Lophophora williamsii or Opuntia ficus-indica (SRR 1616998) unigenes aligned for more than 50 bp with at least 80% blastn similarity. Under these criteria, 50% (14,164) of gene models had evidence of transcription from saguaro RNA-Seq data, and two-thirds (18,597) if including evidence from assembled unigenes.
Robustness of BPP Model Assumptions.
BPP’s substitution model for computing gene tree likelihoods is the Jukes–Cantor (JC) model with a molecular clock. We evaluated model choice in BPP relative to more complex models using likelihood ratio tests (109) with likelihoods computed in PAUP*. Of the 2,291 gene alignments in the most inclusive 50% gene confidence set, the JC model was rejected in favor of HKY85 in 97.8% of the alignments (P = 0.05; df = difference in model parameters = 4; base frequencies estimated). The molecular clock was rejected in far fewer: 18.6–21.6% of alignments (P = 0.05; df = #taxa-2 = 3), depending on the confidence set. Results were similar for other partitions. Although few in number, these genes might conceivably have an outsize effect. To estimate the magnitude of overall lineage effects in Pachycereeae, we concatenated all alignments in the 90% gene set in which the gene tree was concordant with the species tree and estimated branch lengths in PAUP* with an HKY85 + Gamma model (estimating ti/tv ratio). There was evidence of a slight slowdown in rate in the Carnegiea/Pachycereus clade, but overall the effect size for lineage effects was small, with a coefficient of variation in root to tip rates of 13.2%.
To examine model robustness in more detail, we compared population size and divergence time parameter estimates (posterior means) from our data using another program, BEAST 2 v. 2.4.6 (63), which implements more substitution models, including a relaxed clock. We estimated parameters in nine independent sets of 50 genes from the 90% gene All partition set comparing either JC/HKY85 models (assuming a clock), or clock/UCLN (uncorrelated lognormal) models (assuming JC). Input files for BEAST, in XML format, were generated in BEAUti using the StarBeast2 package. Models were unlinked across genes, sites, trees, and clock rates. We estimated a single population size for each branch by setting the population model to “constant” population size. We assigned the same priors in BEAST as we had used in BPP, noting that population sizes in BEAST are in units of 2Neµ, whereas BPP uses units of 4Neµ. “Clock rates” and SDs for each gene tree were unlinked, and we assumed the default prior exponential distribution for the SD. Number of generations was set to 10–50 million as needed to obtain ESS values above 200; information from the first 10% of the MCMC chains was discarded.
Paired t tests indicated no significant differences between JC and HKY85 models for population size, although estimates at one internal node were nearly significant. Effect sizes of the models (the difference between the mean estimate across replicates in the two treatments) were less than 2.8% (Table S8). Divergence time estimates, which were scaled to the root age, were not significantly different between JC/HKY85 models. Population size estimates were more sensitive to clock/UCLN model differences, but the effect sizes were still small. For two of the three internal edges, population size estimates varied by as much as 10.6% between the clock and UCLN models. Divergence time estimates were less sensitive to the clock/UCLN model assumptions, with no significant t tests and effect sizes of less than 3.4% of the root to tip distance in the tree. Overall, these experiments suggest that these datasets are large enough to uncover significant but subtle differences in parameter estimates between some models. However, these effects are unlikely to be large enough to affect our conclusions.
Materials and Methods
Sampling, Genome Sequencing, and Annotation.
For details of taxa sampled, library construction, genome sequencing, assembly, and annotation, see SI Materials and Methods.
Phylogenomic Analyses.
Gene alignments and gene trees.
Sets of gene trees were inferred from gene alignments constructed from the genome assemblies with custom PERL scripts. CDS and intron regions were extracted from the saguaro genome based on its annotation. Each region was used as a query in blastn searches (72) against the other four cactus assemblies, with an e value cutoff of 10−10 for the intron searches. To enrich for orthologs, (i) a CDS region was kept for further processing if it returned exactly one hit for all four taxa and (ii) the gene was kept if and only if at least two CDS regions from the same gene passed test (i).
We then used predicted saguaro amino acid sequences within the retained CDS regions in pairwise tblastn runs against the nucleotide hits found in each of the four cacti. Thereby, we obtained subject sequences in the same frame as the saguaro query. We used Muscle (73) to align these protein regions and then tranalign (74) to align the underlying nucleotide sequences. Each gene thus consists of CDS regions that have sequence data present for all taxa. Given at least one such region, introns were then evaluated. An intron region was kept if it returned exactly four hits, one per taxon, and each hit covered at least 50% of the query saguaro sequence. Any alignment in which at least one taxon had more than 50% gaps or missing data was excluded from further consideration.
Three partitions of each gene nucleotide alignment were prepared: one having just third positions in codons (“Third”), one having third positions plus introns (“Neutral”); and one having all nucleotide positions (“All”). Bootstrap maximum likelihood majority rule trees for the All positions alignment were constructed with PAUP* v. 4.0a with an HKY85 model (75). “Gene confidence sets” were compiled based on the quality of their trees (e.g., the “90%” set is all genes for which the minimum bootstrap value is 90% for all clades). Sets were assembled at 50, 60, 70, 80, and 90% levels and rooted with Pereskia.
Species tree inference.
Species trees were constructed first by an alignment-based Bayesian method, BPP v.3.3a (35, 61), in which the input was the alignments from a gene confidence set and partition. Convergence of Markov chains was checked by running three chains from random starting trees for each of the five gene confidence sets (burnin = 4,000; number of generations = 40,000). This is aimed at avoiding multiple optima in the posterior (76) and may be more informative than examining acceptance rates or effective sample sizes (35).
The prior for the scaled root divergence time, τroot, was set to a gamma distribution, based on the mean and variance of pairwise Ks distances between saguaro and Pereskia (Table S7). The prior for scaled population size, θroot, was set to a gamma distribution with a mean estimated from the modern nucleotide diversity of saguaro, which was estimated, using the program PSMC (77) applied to the saguaro genome scaffolds ≥100 kb in length, to be π = 0.00146 (variance set to half the mean). The rate variation prior was set to a Dirichlet prior with parameter α = 2.
A second, tree-based pseudolikelihood method, MP-EST v. 1.5 (34), was used with an input of all rooted gene trees constructed from the All partition from each gene confidence set. In each, five replicate searches from random starting species trees were done. To check for multiple optima, an exhaustive enumeration of pseudolikelihood scores for all possible ingroup rooted species trees was done by supplying MP-EST with those 15 trees.
Gene tree discordance.
Discordance between gene trees and the inferred species tree was assayed on a clade basis and for the entire gene tree jointly using PAUP*. Discordant gene trees were binned into groups of identical topology by comparing them to all 14 possible rooted discordant binary trees for five taxa using a custom PERL script. Topological discordance in the All partition was visualized using DensiTree (78) by making the gene trees in the 90% gene confidence set ultrametric using a clock model in PAUP* and scaling root ages to 1.0 (Fig. 2A).
Introgression and network reconstruction.
We used InferNetwork_ML in PhyloNet (36) to reconstruct optimal phylogenetic networks using maximum likelihood, using as input gene trees from the All partition [invoked with “InferNetwork_ML (all) h -n 10 -di -o -po -x 20 -s starting_tree”, where h is the maximum allowed number of reticulation events]. Values of h of 0 (a tree) and 1 were each run from three different starting topologies. Each inferred reticulation node has two incident edges with inheritance probabilities γ and 1 − γ. The tree imbedded in the network, obtained by following the path along the larger of the two inheritance probabilities, is called the major tree.
To compare solutions, we used the AIC information criteria (37): AIC = 2k − 2 log L, where k is the number of parameters in the model and L is the maximum likelihood value of the model. To assess the statistical support for the optimal network, we used PhyloNet’s parametric bootstrap procedure with the same search parameters as used previously.
Demographic inference.
We used BPP to infer divergence times and effective population sizes conditional on the species tree found above. To ensure convergence of the Markov chains, we increased the number of generations to 400,000 (burnin = 40,000), leading to all parameter estimates having an acceptable effective sample size (ESS) >500 (79) and acceptance rates lying in the range 0.25–0.40. We also examined parameter traces in Tracer (63). The raw parameter output of BPP is in units of sequence divergence, D = Tμ, where T is the time in generations and μ is the per generation mutation rate per site, and scaled effective population size, θ = 4Neμ. Note that 2D/θ = T/2Ne is the edge length in “coalescent time units.”
Quantifying Hemiplasy.
For any site in an alignment, let mS and mG be the inferred number of state changes on the species and gene trees respectively. If the two trees are concordant, mS = mG, but if the trees are discordant, then the number of changes may differ. The simplest case of hemiplasy is mG = 1 but mS > 1: A homology on the gene tree is a homoplasy on the species tree (Fig. 1, Lower Left). A site may also exhibit homoplasy on the gene tree: mG > 1 (Fig. 1, Lower Right). We define a hemiplastic site as one in which mS > mG. The exceptional case of mG > 1 but mS = mG is “homoplasy but not hemiplasy.”
Ancestral states for each residue in each protein sequence alignment in the most robust 80% and 90% gene confidence sets were reconstructed using parsimony with PAUP*, and mS and mG were computed for all potentially parsimony informative sites. Sites homoplastic on the species tree for concordant gene trees, and sites homoplastic or hemiplastic for nonconcordant gene trees, were tallied by a custom PERL script.
Neutral Mutation Rate Estimates.
We estimated neutral mutation rate from BPP output as sequence divergence from root to tip divided by the crown group age of cacti, since the assumed model is ultrametric. We also used codeml (80) to infer pairwise synonymous divergences in the coding regions across the entire set of gene trees, using the crown age of Cactaceae estimated at 26.88 MYR (81).
Acknowledgments
We thank S. Kumar, T. Hernández-Hernández, B. Rannala, K. Steele, F. Tax, L. Venable, and D. Zwickl for discussion, and the Tumamoc Hill Reserve, Tucson, and the Desert Botanical Garden, Phoenix, for permission to collect material. Funding was provided by the University of Arizona–Universidad Nacional Autónoma de México Consortium for Drylands Research, the Tucson Cactus and Succulent Society, and Arizona State University’s College of Liberal Arts and Sciences and the School of Life Sciences. A.B. received sabbatical support at the University of Arizona from DGAPA-Universidad Nacional Autónoma de México and by Programa de Apoyo a Proyectos de Investigación e Innovación Tecnológica, Universidad Nacional Autónoma de México Grant IN213814. N.K.W. received support from NIH Grant R35GM119816.
Footnotes
- ↵1To whom correspondence should be addressed. Email: sanderm{at}email.arizona.edu.
Author contributions: D.C., A.B., M.M.M., N.K.W., R.A.W., M.F.W., and M.J.S. designed research; D.C., A.B., E.B., J.L.M.C., M.M.M., M.F.W., and M.J.S. performed research; D.C., K.L.C., L.E.E., S.L., T.L.L., M.M.M., M.F.W., and M.J.S. analyzed data; and D.C., A.B., E.B., L.E.E., T.L., M.M.M., N.K.W., M.F.W., and M.J.S. wrote the paper.
The authors declare no conflict of interest.
This article is a PNAS Direct Submission.
Data deposition: The sequences reported in this paper have been deposited in GenBank under BioProject number PRJNA318822; individual accession numbers are listed in Tables S1, S3, and S5.
This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.1073/pnas.1706367114/-/DCSupplemental.
Published under the PNAS license.
References
Genome evolution in columnar cacti
Article Classifications
- Biological Sciences
- Evolution
Sign up for Article Alerts
Jump to section
You May Also be Interested in
Researchers are exploring whether these ubiquitous fluorinated molecules might worsen infections or hamper vaccine effectiveness.
Image credit: Shutterstock/Dmitry Naumov.
One protein seems to play a key role in touch sensitivity for flytraps and other meat-eating plants.
Image credit: Shutterstock/Kuttelvaserova Stuchelova.
Adam Mastroianni and Daniel Gilbert explore why conversations almost never end when people want them to.
A study finds that giant pandas roll in horse manure to increase their cold tolerance.
Image credit: Fuwen Wei.