Gene family evolution and natural selection signatures in Datura spp. (Solanaceae)

Introduction

One of the fundamental endeavors in biology is to unravel the process and mechanisms of species diversification (Savolainen et al., 2013). The development of phylogenetic theory (Hennig, 1950, 1966), followed by the generation of massive DNA sequence datasets, increased computing power, and the proliferation of analytical methods (e.g., maximum likelihood, Bayesian inference, Bayesian molecular dating), has led to the recognition of the temporal dynamics of evolutionary radiation (Magallon and Sanderson, 2001; Alfaro et al., 2009; Stadler, 2011; Nürk et al., 2020). In addition, recent advances in genomics and bioinformatics have allowed the elucidation of the transcriptome and genome assemblies of several plant species (Leebens-Mack et al., 2019), which constitute the basic resources to find candidate genes underlying physiological, developmental, and ecological traits involved in species diversification (Eckardt, 2001; Ehrenreich and Purugganan, 2006; Anderson et al., 2011; Chang et al., 2012; Savolainen et al., 2013; Flood and Hancock, 2017; Sork, 2018; Bamba et al., 2019; Kariñho-Betancourt et al., 2020; Rich-Griffin et al., 2020). Genomic approaches that unravel physicochemical divergence, natural selection, and/or gene expansions/contractions in gene families of related taxa have also contributed to enhanced understanding of the evolutionary trajectories of specific traits (Stone and Sidow, 2005; Arendt and Reznick, 2008; Han et al., 2013; Murrell et al., 2013; Emms and Kelly, 2019). For instance, it is possible to discover candidate genes involved in the production of secondary plant metabolites that could also play important roles in mediating plant-based interactions in the evolution of differentiated defensive mechanisms (Bergelson et al., 2001; Friedman and Baker, 2007; Ramakrishna and Ravishankar, 2011; Chang et al., 2012; Kariñho-Betancourt, 2018; Kourelis and van der Hoorn, 2018; Wink, 2018; Erb and Kliebenstein, 2020; Rich-Griffin et al., 2020).

The comparative genomics study (e.g., gene family analysis) of plants has provided a broad view of the mechanisms at the molecular level behind species evolution and diversification (e.g., genomic signatures of directional selection on particular genes) (Song and Mitchell-Olds, 2011; Savolainen et al., 2013; Flood and Hancock, 2017; Russell et al., 2017; Anderson and Song, 2020). Genomic studies have been central to investigating the evolution and diversification of ecologically and economically relevant plant species, such as the Solanaceae family (Soltis and Soltis, 2021). Several species of this plant family have been subjected to metabolomic and genomic studies on specialized secondary metabolites, such as terpenoids, tropane alkaloids, glycoalkaloids, and phenolics, implicated as defenses against natural enemies (Friedman, 2002, 2006; Cheng et al., 2007; Sierro et al., 2014; Akyol et al., 2016; Chowański et al., 2016; Xu et al., 2017; De-la-Cruz et al., 2021). However, this is not the case for nonmodel plant species in the Solanaceae family – such as Datura species – where few genomic studies examining the processes and candidate genes that underlie evolutionary diversification have been carried out.

Within the Solanaceae family, the species of the genus Datura produce high concentrations of tropane alkaloids (Jirschitzka et al., 2012; Kohnen-Johannsen and Kayser, 2019; Céspedes-Méndez et al., 2021), metabolites long used as medicines (Roddick, 1991; Geeta and Gharaibeh, 2007; Benítez et al., 2018; Sahu et al., 2022). Ecological evidence indicates that some of these chemical compounds, as well as some terpenoids, play central roles as defenses against herbivores and possibly against pathogens, viruses, and fungi in Datura (Shonle and Bergelson, 2000; Hare and Elle, 2001; Hare and Smith, 2005; Castillo et al., 2014; Kariñho-Betancourt et al., 2015; Miranda-Pérez et al., 2016; De-la-Cruz et al., 2020a,2021; Goldberg et al., 2020; Céspedes-Méndez et al., 2021; Velázquez-Márquez et al., 2021). Furthermore, there is evidence that these species produce different concentrations of total tropane alkaloids (Doncheva et al., 2006; Kariñho-Betancourt et al., 2015). The genus Datura includes 13–14 species of annual herbs and perennial shrubs native to dry, temperate, and subtropical regions of North America, distributed mainly in Mexico, which is considered its center of origin (Avery, 1959; Barclay, 1959; Symon and Haegi, 1991; Luna-Cavazos and Bye, 2011). Plant size ranges from 30 cm to 2 m in height (Avery, 1959). Flowers open at dusk and are pollinated by hawkmoths (Jiménez-Lobato and Núñez-Farfán, 2021). All species have tubular erect flowers with large variations in size and color (Matuda, 1952). While most ruderal species have small-sized white, yellow and purple flowers (3–8 cm), plants distributed in less disturbed environments have mostly large white flowers (10–17 cm). Nevertheless, the genomic mechanisms underlying the evolutionary divergence of Datura species are unknown.

In this study, we selected four Datura species of two different clades of the genus (Mace et al., 1999; Bye and Sosa, 2013; Kariñho-Betancourt et al., 2015), some of which also overlap with respect to their distribution ranges within Mexico (Supplementary Figure 1). Datura inoxia, D. wrightii, and D. pruinosa are nested within the section of nodding fruits and larger flowers (except for D. pruinosa). Datura inoxia occurs in central and northern Mexico, D. wrightii occurs in northern Mexico, and Datura pruinosa is distributed in southern Mexico. In contrast, D. stramonium is nested within the section of erect fruits and is distributed throughout North and South Mexico (Mace et al., 1999; Bye and Sosa, 2013; Kariñho-Betancourt et al., 2015). All Datura species studied here are summer annuals, except for D. wrightii, which is perennial. Hence, these species represent good systems for addressing their evolutionary divergence at the genomic level. To this end, we performed a thorough gene family analysis using de novo leaf transcriptomes of these four Datura species to find candidate genes evolving under positive selection, physicochemical divergence (changes in physicochemical properties of proteins) and expansion (gain of genes).

Materials and methods

Plant material and sequencing

Our sample design consisted of nine plants per species (different maternal families) from the four Datura species (see above). Seeds were collected in the field (natural progenies; Supplementary Figure 1) and germinated in the glasshouse of the Institute of Ecology, UNAM (UTM 19.3218, −99.1922). Once true leaves appeared and seedlings reached 3–5 cm long, they were planted in 150 mL pots in sterilized soil, watered ad libitum, and grown under a 16:8 L:D cycle at 25:20°C (L:D). Then, we sampled fully expanded leaves from each plant of the four Datura species before and during flowering (dx.doi.org/10.17504/protocols.io.bx4zpqx6). All leaves were flash-frozen in liquid nitrogen and kept in a −80°C freezer. We extracted total RNA from the leaves of each plant species of Datura studied here using the TRIzol extraction method (doi.org/10.17504/protocols.io.bx4zpqx6) (Rio et al., 2010). RNA quality and quantity were determined using a Nanodrop 2000 instrument (Thermo Scientific) and Bioanalyzer Chip RNA 7500 series II (Agilent). Each sample was sequenced using an Illumina NextSeq 500 in paired-end mode (2 × 75 bp). The library preparation and sequencing were carried out at the Unidad Universitaria de Secuenciación Masiva y Bioinformática (UUSMB) of the Institute of Biotechnology, UNAM. The quality of the RNA sequences of each plant individual was verified using FastQC (Andrews, 2010). Trimmomatic (Bolger A. M. et al., 2014) was used to remove sequences below 20 Phred quality scores.

De novo transcriptome assembly and validation

We pooled the sequences of all analyzed plants of each Datura species as input for a de novo assembly using Trinity v2.12.0 (Grabherr et al., 2011). After transcriptome assembly, we retained only the longest isoforms using TransDecoder v5.5.0.¹ A further step was carried out using CD-HIT v4.8.1 to cluster similar proteins and remove their redundancy (Li and Godzik, 2006; Fu et al., 2012) to retain only “primary transcripts” (Li and Godzik, 2006; Fu et al., 2012). We evaluated the transcriptome assemblies using the standard metrics (total genes and transcripts, percent of GC, Nx length statistics and median and average contig) with the script “TrinityStats.pl” of Trinity v2.11.0 (Grabherr et al., 2011). We assessed transcriptome completeness using the Benchmarking Universal Single-Copy Orthologs (BUSCO), which determines completeness by determining the fraction of a lineage-specific set of conserved genes present in an assembly. We performed the analysis using BUSCO v.5.1.2 (Seppey et al., 2019) with the transcriptome mode option. The database “embryophyta_odb9” was used as a reference for searching BUSCOs in the four Datura transcriptomes (Simão et al., 2015). BUSCOs were classified as complete and single-copy (S), complete and duplicated (D), fragmented (F), or missing (M).

Ortholog identification and phylogenetic inference

Four Datura and three Nicotiana transcriptomes were used for phylogenetic and gene family evolution analyses. Since we encountered logistic issues downloading other Solanaceae transcriptomes (e.g., Solanum transcriptomes), we decided to use these three Nicotiana species for gene family construction (see De-la-Cruz et al., 2021). Thus, protein-coding genes and CDS from the three Nicotiana transcriptomes were sourced from the Sol Genomics Network² [N. tomentosiformis (Sierro et al., 2013), N. sylvestris (Sierro et al., 2013), and N. attenuata (Xu et al., 2017)]. We used the combined set of Datura and Nicotiana proteomes to infer orthogroups in the OrthoFinder v2.5.2 program (Emms and Kelly, 2019). Briefly, DIAMOND blast (Buchfink et al., 2015) was used in OrthoFinder (Emms and Kelly, 2019) for orthogroup inference (E-value < 1e^–5) and the MCL algorithm was employed for clustering of similar sequences (Enright et al., 2002). MAFFT v7 was used as a multiple protein sequence aligner of each orthogroup/gene family (Katoh et al., 2002), and FastTree2 v2.1.10 (Price et al., 2010) was used for maximum likelihood gene tree inference. The species tree inference was constructed using a concatenated alignment of 3,219 single-copy orthogroups in OrthoFinder (Emms and Kelly, 2019). The species tree was inferred with FastTree2 (Price et al., 2010). Rooting was performed using the STRIDE algorithm (Species Tree Root Inference from Duplication Events; Emms and Kelly, 2019) of OrthoFinder, which selected N. sylvestris, N. attenuata, and N. tomentosiformis as outgroups of the Datura phylogeny.

Ultrametric tree and gene family evolution analysis

To build an ultrametric phylogeny (i.e., estimating divergence times and rates of substitution in a phylogeny) for the analysis of gene family evolution (expansions and contractions in gene families; see below), we used the rooted species tree obtained from OrthoFinder to search the divergence times between the branches using the TimeTree webtool (Kumar et al., 2017). Then, the rooted species tree and the information of divergence times were used to create an input file using the script ‘‘cafetutorial_prep_r8s.py’’ from the CAFE tutorial webpage³ to run the r8 v1.5 program for ultrametric inference (Sanderson, 2003). The tip-to-root length was adjusted to match the approximately 35-million-year evolutionary history between Nicotiana and Datura species (Bombarely et al., 2016; Kumar et al., 2017; De-la-Cruz et al., 2021).

To evaluate the rate of gene gain and loss across the Datura and Nicotiana species studied here, we used only the gene families with more than three genes per family (24,884) and the ultrametric species tree as inputs to CAFE v4.2.1 (Han et al., 2013). CAFE (Computational Analysis of gene Family Evolution) estimates the birth and death parameter (λ) for the provided tree and gene family counts. The λ parameter describes the probability that any gene will be gained or lost (Han et al., 2013). First, the Python scripts provided by the CAFE pipeline were used to estimate the error in our dataset (Han et al., 2013). Then, CAFE was run using the mode in which the gain and loss rates were estimated together (λ) for the whole phylogeny (Han et al., 2013). The CAFE overall p value threshold was kept at its default value (0.01) for the entire analysis. We used the script “cafetutorial_report_analysis.py” from CAFE to obtain summary statistics. The script parseCafeResult.R from the SlydGeneFamsAnalyses R package (see De-la-Cruz et al., 2021) was used to parse the CAFE outputs for functional enrichment analysis (see below).

Multivariate analysis of protein polymorphisms

Multiple sequence alignments of the 24,884 protein families (families with more than three proteins) were used to carry out a multivariate analysis of protein polymorphisms (MAPP program) (Stone and Sidow, 2005). MAPP quantifies the physicochemical variation (hydropathy, polarity, charge, volume, free energy in alpha-helix conformation, and free energy in beta-strand conformation) in each column of a multiple sequence alignment and calculates the deviation of candidate amino acid replacements from this variation. The greater the deviation, the higher the probability that a replacement impairs the function of the protein and the greater its predicted effect on the protein’s function (Stone and Sidow, 2005). Therefore, MAPP was used to estimate the physicochemical divergence in each protein family. First, we used the script readAndParseOrthogroupsTxt.R from the SlydGeneFamsAnalyses R package (De-la-Cruz et al., 2021) to parse and create folders from each protein family and store the corresponding protein tree and multiple sequence alignment from OrthoFinder results. Then, the MAPP program was run with default parameters for each of the protein families (Stone and Sidow, 2005). The script readMappResults.R in the SlydGeneFamsAnalyses R package (De-la-Cruz et al., 2021) was used to parse and read all the MAPP results of gene families. This script reads the MAPP results for all gene families, adjusts the p value, finds genes of families with good multiple sequence alignments (Valdar Score > 0.6), and only retains the significant sites with physicochemical divergence that correspond with conserved domain proteins (De-la-Cruz et al., 2021). This method allows us to score residues in a multiple sequence alignment and assigns a score ranging from 0 for low conservation to 1 for high conservation (Valdar, 2002).

Divergent positive selection in gene families

A codon-level analysis of positive natural selection was performed with FUBAR (Fast, Unconstrained Bayesian Approximation) (Murrell et al., 2013) on 24,884 gene families. FUBAR uses a Bayesian approach to infer nonsynonymous (dN) and synonymous (dS) substitution rates on a per-site basis for a given coding alignment and corresponding phylogeny (Murrell et al., 2013). This method assumes that the selection pressure on each site is constant along the entire phylogeny (Murrell et al., 2013). FUBAR reports evidence for positive selection using posterior probabilities (which range from 0 to 1), not p values. Generally, posterior probabilities > 0.9 strongly suggest positive selection (Murrell et al., 2013).

To run FUBAR, we first removed trailing stop codons from the CDS in each transcriptome and then applied PAL2NAL (Suyama et al., 2006) to produce a codon alignment for each gene family. PAL2NAL is a program that converts a multiple sequence alignment of proteins and the corresponding DNA (CDS) sequences into a codon alignment (Suyama et al., 2006). Therefore, we used the protein tree that we already had from each protein family to run PAL2NAL. FUBAR was run for each protein family’s codon alignments (De-la-Cruz et al., 2021). A custom Python script was used to transform the “.json” format from the FUBAR result to tabular format. Then, the R script “loadFubarResults.R” from the R package GeneFamilies (De-la-Cruz et al., 2021) was used to obtain a table with the significant posterior probabilities of a codon being subject to positive selection for each gene family (significant posterior probabilities ≥ 0.98; Bayes Factor > 100). As stated above, we only used alignments with a Valdar score > 0.6.

Enrichment analysis

Fisher’s exact test (Fisher and Russell, 1922) was employed for functional enrichment/overrepresentation analyses to look for biological functions enriched in a gene list to a greater extent than expected by chance. In other words, the overrepresentation method compares the set of all genes annotated to pathways to a list of those genes that are significantly expanded, positively selected and which exhibit physicochemical divergence between species. To this end, we used all the proteins from the seven transcriptomes (262,248 proteins) as the background to detect overrepresented proteins that showed signals of expansion, physicochemical divergence (MAPP), and positively selected conserved amino acids (codons) (FUBAR) (De-la-Cruz et al., 2021). Functional annotation of all transcriptomes was performed using MapMan4 (Schwacke et al., 2019). The enrichment analysis was carried out using our custom script “enrichedAnnosInExpContrFams.R (CAFE)”, “identifyDomainsAtSelectedSites.R (FUBAR)” and “readMappResults.R (MAPP)” from the SlydGeneFamsAnalyses R package (De-la-Cruz et al., 2021). In addition, protein domains were annotated for the four complete Datura transcriptomes using Interproscan v.5.24 (Jones et al., 2014).

Results

Transcriptome sequencing and assembly

Illumina sequencing yielded between 1 and 14 M (million) paired-end reads per sequenced individual of every Datura species (Table 1 and Supplementary Table 2). This resulted in good-quality RNA-seq datasets containing between 86 and 104 M paired-end reads for each species (Table 1 and Supplementary Table 1). The raw read data from Illumina sequencing for each species are deposited in the NCBI under BioProject PRJNA669339. Trinity constructed 81,748 transcripts for D. inoxia, 103,954 for D. pruinosa, 113,343 for D. stramonium and 111,196 for D. wrightii (Table 1). After polishing with TransDecoder and CD-HIT, we retained only 26,112, 30,191, 33,851, and 30,385 transcripts for D. inoxia, D. pruinosa, D. stramonium, and D. wrightii, respectively, for subsequent analyses (see below and Table 1). The percentages of GC and N50 were relatively similar among the four Datura transcriptomes (Table 1). The Datura transcriptomes covered 86.1, 91.8, 90.3, and 92.8% complete single-copy BUSCOs for D. inoxia, D. pruinosa, D. stramonium, and D. wrightii, respectively (Supplementary Table 2). Overall, the four assemblies were relatively similar in gene completeness, with D. wrightii and D. inoxia having the highest and lowest BUSCO scores, respectively (Supplementary Table 2).

Table 1. RNA reads produced by the transcriptome assembly of four Datura species.

Protein family construction

Protein coding genes from three Nicotiana species sourced from the Sol Genomics Network (see “Materials and methods” section) and proteomes/CDS of four Datura species were used to construct gene families/orthogroups. From 262,248 genes, OrthoFinder assigned 85.6% to 30,065 orthogroups/protein families. Fifty percent of all genes were in orthogroups with eight or more genes (G50 = 8; Supplementary Table 3). There were 3,219 orthogroups that consisted entirely of single-copy genes and 11,786 orthogroups with all species present (Supplementary Table 3). The two species sharing the most protein families were N. tomentosiformis and N. sylvestris. Likewise, compared to the other Nicotiana species, N. attenuata shares the highest number of orthogroups with the four Datura species (Supplementary Figure 2). Within the Datura genus, D. wrightii and D. pruinosa shared the highest number of orthogroups (protein families), followed by D. stramonium and D. pruinosa. Datura inoxia was the species that shared the lowest number of orthogroups with the other three Datura species (Supplementary Figure 2). This difference may be explained by the lower number of transcripts/genes assembled for D. inoxia.

According to the estimated calibration time, Nicotiana and Datura diverged ∼25 Mya (Figure 1). Our evidence indicated that D. inoxia, D. pruinosa, and D. wrightii diverged from D. stramonium ∼2.6 Mya. Datura inoxia and D. wrightii were the most recently derived species within the Datura group (Figure 1). The rate of gene gain and loss (λ) resulting from CAFE analysis was 0.039 per year for the whole tree (Figure 1). The internal branch with the largest number of protein gains (expansions) corresponded to the most recent common ancestor of N. tomentosiformis (Figure 1). The internal branch with the largest contractions corresponded to the most recent common ancestor of the Datura clade (Figure 1). Within the Datura group, there was no evidence of expansions and contractions for the most recent common ancestor of D. pruinosa, D. inoxia, and D. wrightii. The terminal branch with the largest number of expansions corresponded to N. sylvestris (Figure 1 and Table 2). Within the Datura clade, the terminal branch with the largest number of expansions was the branch leading to D. stramonium. The terminal branch with the greatest number of contractions was the branch leading to N. attenuata (Figure 1 and Table 2). Within the Datura group, the terminal branch with the largest contractions was the branch leading to D. inoxia, followed by D. stramonium (Figure 1 and Table 2). The internal branch with the largest rapidly evolving gene families corresponded to the most recent common ancestor of N. attenuata and N. sylvestris (Figure 1). The terminal branch with the most rapidly evolving gene families was the branch leading to N. attenuata (Figure 1 and Table 2). Within the Datura clade, the terminal branch with the most rapidly evolving gene families was the branch leading to D. inoxia, followed by D. wrightii (Figure 1 and Table 2). According to CAFE, the number of protein families that experienced neither expansion nor contraction in the genomes of all species tested here varies from 15,896 to 19,743 (Table 2).

Figure 1. Phylogenetic time tree constructed from concatenated alignment of 3,219 orthogroups consistent entirely of single-copy genes from four Datura species and three Nicotiana species. Nicotiana and Datura species diverged approximately 25 Mya. Expansions (red), contractions (blue) and rapidly evolving gene families (green) are shown at nodes and terminal branches. The rate of gene gain and lost (λ) resulted from CAFE analysis was 0.039 for the whole tree.

Table 2. Comparative computational analysis of gene family evolution (CAFE analysis) for seven Solanaceae transcriptomes.

Enrichment tests overlapping families with signals of positive selection, expansion, and physicochemical divergence

We found a total of 11, 14, 10, and 11 enriched protein families (i.e., functionally overrepresented protein families) with a signal of expansion in D. inoxia, D. pruinosa, D. stramonium, and D. wrightii, respectively (Figure 2 and Supplementary Table 4). We also detected 22, 17, 12, and 24 enriched proteins with signals of physicochemical divergence (Figure 3 and Supplementary Table 5) and 51, 59, 51, and 54 enriched proteins with positively selected conserved amino acids (Figure 4 and Supplementary Table 6) for D. inoxia, D. pruinosa, D. stramonium, and D. wrightii, respectively. Several of these enriched family proteins with signals of expansion, positive selection, or physicochemical divergence have been related to responses to specific abiotic and biotic pressures, as discussed below. We found that one and 29 expanded and positively selected protein families were shared between all Datura species, while there was no protein family with significant physicochemical divergence shared between all Datura species (Figures 2–4).

Figure 2. The UpSet plot shows the intersections of the set of a particular expanded protein family annotated with Mapman4 from the seven transcriptomes. Filled cells indicate that the set is participating in the intersection and that the species share that protein family with signal of expansion (A). Heatmap of the enriched expanded proteins of the four Datura species studied here (B). The brightness of color indicates the number (counts) of enriched protein families with signal of expansion that are shared between species. See also Supplementary Tables 4–6.

Figure 3. UpSet plot shows the intersections of the set of a particular protein family annotated Mapman4 with signal of physicochemical divergence from the seven transcriptomes. Filled cells indicate that the set is participating in the intersection and that the species share that protein family with signal of physicochemical divergence (A). Heatmap of the enriched proteins with signal of physicochemical divergence of the four Datura species studied here (B). The brightness of color indicates the number (counts) of enriched protein families with signal of physicochemical divergence that are shared between species. See also Supplementary Table 5.

Figure 4. UpSet plot shows the intersections of the set of a particular protein family annotated with Mapman4 with signal of positive selection from the seven transcriptomes. Filled cells indicate that the set is participating in the intersection and that the species share that protein family with signal of positive selection (A). Heatmap of the enriched proteins with signal of positive selection of the four Datura species studied here (B). The brightness of color indicates the number (counts) of enriched protein families with signal of positive selection that are shared between species. See also Supplementary Table 6.

Candidate genes related to chemical defense and plant resistance (R genes)

We detected several protein families related to plant resistance to pathogens, herbivores and viruses. For instance, the UDP-dependent glycosyl transferase protein family - related to cytokinin conjugation - was expanded and positively selected only in D. stramonium (Figure 2). We also detected enriched protein families with signals of positive selection related to glycosyltransferases and oxidoreductases in all Datura species. Expansion of oxidoreductases was also detected only in D. pruinosa and D. stramonium, while expansion of glycosyltransferases occurred in all the species except in D. stramonium. Moreover, we also observed positive selection for terpenoid proteins in D. pruinosa and D. wrightii, but not in D. inoxia and D. stramonium (Figure 4 and Supplementary Table 6). However, we also found physicochemical divergence with respect to 2-C-methyl-D-erythritol 4-phosphate (of the terpenoid family) only in D. stramonium (Figure 3 and Supplementary Table 5).

Genes related to plant immunity, such as phosphorylation-related proteins belonging to the large family of leucine-rich repeat receptor-like protein kinases (LRRs), showed expansion, physicochemical divergence, and positive selection in all Datura species (Supplementary Tables 4–6). We also found several enriched proteins belonging to the kinase superfamily with physicochemical divergence, expansion and positive selection (Supplementary Tables 4–6), such as the SD-1 kinase protein, which was positively selected in all Datura species (Figure 4). Plant nucleotide binding proteins, such as leucine-rich repeat (NLR) receptors, were only positively selected in D. pruinosa and D. wrightii (Figure 4), and proteins related to receptor-like protein kinase (RLCK-IXb) were expanded in D. inoxia (Figure 2 and Supplementary Table 4).

Candidate genes related to abiotic stressors

We also found different enriched proteins with signals of expansion, physicochemical divergence and positive selection in the Datura species studied here (Supplementary Tables 4–6) associated with abiotic stresses. For instance, we found physicochemical divergence in a protein family related to glutamate biosynthesis (the pyrroline-5-carboxylate synthase gene family) in D. pruinosa (Figure 3). In addition, we found physicochemical divergence in different protein families related to solute transport, such as the metal chelator transporter TCR and metal chelator nicotianamine transporters in D. pruinosa (Figure 4). Datura stramonium showed expansion in the class of Tau glutathione S-transferases (Figure 2 and Supplementary Table 4). In the case of D. inoxia, we detected a signal of expansion in receptor-like cytoplasmic kinases (RLCKs) (Figure 2 and Supplementary Table 4), while we did not find expansions in TCR, RLCKs, and Tau glutathione S-transferases in D. wrightii.

Discussion

In this study, the completeness of single-copy orthologs and general statistics (N50, number of genes) suggest a good quality assembly for all de novo Datura transcriptomes (Simão et al., 2015). The number of complete BUSCOs for all transcriptome assemblies is very similar to that previously reported for D. stramonium and other Solanaceae species, such as tomato and potato, Nicotiana and its wild relatives, as well as Petunia inflata and Petunia axilaris (Xu et al., 2011, 2017; Sato et al., 2012; Sierro et al., 2013; Bolger A. et al., 2014; Bombarely et al., 2016; Hulse-Kemp et al., 2018; Razali et al., 2018; Barchi et al., 2019; De-la-Cruz et al., 2021). Our comparative gene family analysis also revealed signals of positive selection, physicochemical divergence and expansions for candidate genes possibly associated with species divergence.

Our results revealed that the divergence times reported here are consistent with other Solanaceae phylogenies (Särkinen et al., 2013; Bombarely et al., 2016; De-la-Cruz et al., 2021). We also found that N. attenuata shares the highest number of orthogroups with the four Datura species. This result could suggest a closer evolutionary history (i.e., genes forming orthogroups due to a shared origin from a single ancestral gene of the last common ancestor of the compared genomes) between these species (Koonin, 2005). Nevertheless, duplication events could also generate many of these genes forming orthogroups (Koonin, 2005). We detected positive selection of the same 29 protein families of all Datura species. This result could suggest parallel evolution of these genes among all Datura species.

As in many other plant species, we also found several enriched protein families (significantly, in more from the full protein set) in the Datura transcriptomes, with signals of expansion, physicochemical divergence and/or positive selection varying among species. For instance, UDP-dependent glycosyl transferase was expanded and positively selected only in D. stramonium. In plants, UDP-glycosyltransferases glycosylate various phytohormones and metabolites (e.g., triterpenoids) in response to biotic and abiotic stresses such as defenses to herbivores (Vogt and Jones, 2000; Rehman et al., 2018; Rahimi et al., 2019; Wilson and Tian, 2019). Indeed, we also detected enriched protein families with signals of positive selection related to glycosyltransferases and oxidoreductases in all Datura species. Transferases and oxidoreductases are enzymes directly involved in the biosynthesis of terpenes, tropane alkaloids (Kanehisa and Goto, 2000; Dräger, 2006; Jirschitzka et al., 2012; Kanehisa et al., 2016, 2017) and many other plant defensive compounds, such as phenolics, glucosinolates, salicylates, and anthocyanins, via glycosylation (Vogt and Jones, 2000; Bowles et al., 2005; Wilson and Tian, 2019). Expansion of oxidoreductases was also detected only in D. pruinosa and D. stramonium, while expansion of glycosyltransferases occurred in all the species except D. stramonium. Moreover, we also found positive selection of terpenoid proteins in D. pruinosa and D. wrightii, but not in D. inoxia and D. stramonium (Figure 4 and Supplementary Table 6). However, we also found physicochemical divergence in the protein family related to 2-C-methyl-D-erythritol 4-phosphate (terpenoids) only in D. stramonium (Figure 3 and Supplementary Table 5). This molecule is one of the two major building blocks of terpenoid compounds (Rohdich et al., 2000; Narayanasamy et al., 2010). Thus, it is likely that this divergence in the physical and chemical properties in this protein family could be related to changes in the production and/or concentration of specific terpenoid compounds. A recent study carried out in natural conditions with experimental plants of D. stramonium found that at least one triterpenoid compound was positively selected (Darwinian fitness) as a defense against the larvae of the specialized beetle Lema trilienata daturaphila, which is responsible for most of the leaf damage to plants (De-la-Cruz et al., 2020b). Thus, we speculate that terpenoid genes (and their products) could be involved in the production of defensive compounds to guard against natural enemies (e.g., as we have found in D. stramonium).

Terpenoids are involved in several biological functions related to development and defense against herbivores, fungi and pathogens (Degenhardt et al., 2009; Singh and Sharma, 2015; Boncan et al., 2020; Zhou and Pichersky, 2020). Moreover, several terpenoids are volatilized to attract predators and parasitoids of herbivores (Singh and Sharma, 2015; Lackus et al., 2018; Boncan et al., 2020). Since the four Datura species are attacked by the specialist herbivore L. trilineata daturaphila (Kogan and Goeden, 1970; Castillo et al., 2014; De-la-Cruz et al., 2020b; Goldberg et al., 2020), we speculate that the positive selection (D. pruinosa and D. wrightii), expansion (only in D. wrightii) and physicochemical divergence (only in D. stramonium) detected for terpenoids could be related how these species cope with their natural enemies such as herbivores, pathogens and viruses (Pichersky and Raguso, 2018; Karunanithi and Zerbe, 2019; Boncan et al., 2020; De-la-Cruz et al., 2020a,b; Goldberg et al., 2020, 2021; Zhou and Pichersky, 2020; Velázquez-Márquez et al., 2021). Thus, positive selection of terpenoid proteins in almost all Datura species could suggest species-specific evolutionary dynamics related to the different biotic and abiotic pressures faced by each Datura species in their natural environments (Castillo et al., 2014; Kariñho-Betancourt et al., 2015; Pichersky and Raguso, 2018; Karunanithi and Zerbe, 2019; De-la-Cruz et al., 2020b; Goldberg et al., 2020, 2021; De-la-Cruz et al., 2021).

We also found that the large family of leucine-rich repeat receptor-like protein kinases (LRRs) showed expansion, physicochemical divergence, and positive selection in all Datura species. LRR genes play a central role in plant innate immunity and are involved in the regulation of jasmonate/ethylene and salicylate hormones implicated in the defensive response against pathogens and herbivores (McHale et al., 2006; Padmanabhan et al., 2009; Magalhães et al., 2016; Tang et al., 2017). Notably, D. stramonium had the highest number of protein families related to phosphorylation under positive selection. We also found several enriched proteins belonging to the kinase superfamily with physicochemical divergence, expansion and positive selection, such as the SD-1 kinase protein, which was positively selected in all Datura species. Kinase superfamily proteins have been related to plant defense signaling in response to different stresses, including pathogen invasion (Stone and Walker, 1995; Afzal et al., 2008). Likewise, plant nucleotide binding and leucine-rich repeat (NLR) receptors are linked to defensive responses against pathogens triggering plant immunity effectors (Cui et al., 2015; Thordal-Christensen, 2020; van Wersch et al., 2020). In our study, NLR proteins were only positively selected in D. pruinosa and D. wrightii, and proteins related to the receptor-like protein kinase (RLCK-IXb) were expanded in D. inoxia.

On the other hand, we observed a signal of physicochemical divergence in genes related to transcriptional regulation (transcriptional corepressor LUG) among the three Datura species. It has been reported that LUG plays an important role in leaf blade outgrowth and flower development, also enhancing the polarity and growth defects in leaves of Arabidopsis thaliana (Stahle et al., 2009; Zhang et al., 2019). It could be possible that the physicochemical divergence detected in LUG could be involved in the mechanisms determining leaf shape and flower size of these Datura species. Interestingly, the transcriptional corepressor LUG is not enriched in D. pruinosa, which has smaller leaves and flowers than the four Datura species studied here (Figure 4; Matuda, 1952). Thus, LUG genes could be candidates for studying leaf and flower shape and size in Datura.

We detected physicochemical divergence in protein families related to hemicellulose function, such as xylosyltransferase IRX9 genes, among D. inoxia, D. stramonium, and D. wrightii, but not for D. pruinosa. It has been reported that mutations in IRX9 cause a collapsed xylem phenotype and decreases in xylose and cellulose in cell walls, suggesting that IRX9 is required for normal glucuronoxylan elongation (Peña et al., 2007; Wu et al., 2010). Glucuronoxylans occur in angiosperm secondary walls and are the main hemicellulosic polysaccharide in the secondary cell walls of plants, and it has been detected that IRX9 is specifically expressed in cells undergoing secondary wall thickening (Peña et al., 2007; Wu et al., 2010). Moreover, a study in A. thaliana revealed that IRX mutants had reduced cellulose content in the stems, resulting in decreased stiffness of stem material (Turner and Somerville, 1997). Thus, the IRX9 protein is also a candidate gene to study in terms of involvement in leaf thickness, stem stiffness, stem architecture and elongation differentiation in the Datura species studied here.

We also found different enriched proteins with signals of expansion, physicochemical divergence and positive selection in D. pruinosa that have been associated with abiotic stresses. For instance, we observed physicochemical divergence in a protein family related to glutamate biosynthesis (the pyrroline-5-carboxylate synthase gene family) in D. pruinosa. In Datura stramonium, De-la-Cruz et al. (2020b) found that pyrroline, a derivative compound from pyrroline-5-carboxylate synthase, is involved in defense against the folivore beetle Epitrix parvula (Chrysomelidae). This gene is also involved in proline accumulation, a widespread metabolic response of plants to osmotic stress (Pérez-Arellano et al., 2010). In plants, proline is useful not only during conditions of drought and high salt osmotic stress, but also as a response to environmental changes such as low temperature, nutritional deficit, heavy metals, UV radiation, and bacterial pathogens (Parvanova et al., 2004; Haudecoeur et al., 2009; Pérez-Arellano et al., 2010; Qamar et al., 2015). In some species, proline also accumulates during the development of flowers, pollen, ovules or fruits (Venekamp and Koot, 1988; Fujita et al., 1998; Pérez-Arellano et al., 2010).

In addition, we found physicochemical divergence in different protein families related to solute transport, such as the metal chelator transporter TCR and metal chelator nicotianamine transporters in D. pruinosa. Transition metals (e.g., iron, copper, and zinc) are essential for many physiological processes but can be toxic at elevated levels (Clemens, 2001). TCR proteins contribute to metal detoxification by buffering cytosolic metal concentrations (Clemens, 2001). In the same way, D. stramonium showed expansion in the Tau glutathione S-transferase class. The Tau and Phi classes of transferases are the largest groups in plants and play crucial roles in the remediation of environmental pollution by organic xenobiotics, including herbicides and industrial chemicals, as well as salinity, oxidative stress, and drought (Kilili et al., 2004; Benekos et al., 2010; Dixon and Edwards, 2010; Sharma et al., 2014; Gallé et al., 2019). In the case of D. inoxia, we detected a signal of expansion in receptor-like cytoplasmic kinases (RLCKs), which regulate plant cellular activities in response to biotic/abiotic stresses and endogenous extracellular signaling molecules (Liang and Zhou, 2018). Datura plants grow in highly human-disturbed, contaminated habitats and degraded soils, including landfills.

Interestingly, TCR, RLCKs, and Tau glutathione S-transferases are not expanded in D. wrightii, which occurs in less disturbed habitats. Thus, it could be possible that TCR, RLCKs and Tau glutathione S-transferases could be involved in the mechanisms related to the assimilation, long-distance transportation and detoxification of metal compounds in D. pruinosa, D. inoxia, and D. stramonium, respectively, to survive and grow in highly disturbed and polluted habitats. All of these hypotheses should be tested through functional analysis in further studies.

Conclusion

The transcriptome assemblies of these four Datura species provide a valuable molecular resource, as there are few sequenced transcriptomes of nonmodel species from the Solanaceae family. Our results provide new insights into the diversification among these four Datura species. Understanding the evolution of genes involved in the synthesis of secondary compounds and genes related to defense and other functions, such as development and growth, is essential to elucidating their roles in species divergence. Here, we described how the Datura gene families expanded and revealed positive selection and physicochemical divergence of proteins related to terpenoids, R genes, tropanes and other secondary metabolites putatively relevant to the diversification processes among these species. Gene family analysis is a useful tool for understanding the genomic evolution of economically and ecologically significant taxa. It is important to note that the overrepresented candidate genes with signals of positive selection and physicochemical divergence should be tested with respect to functional analysis to evaluate their ecological and evolutionary significance.

Data availability statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi.nlm.nih.gov/, PRJNA669339.

Author contributions

ID-L-C and EK-B conceived and designed the study and wrote the first draft of the manuscript. EK-B designed and performed the experiments, RNA extraction, and data compilation. ID-L-C analyzed the data. KO, JN-F, EK-B, and ID-L-C contributed reagents, materials, and analysis tools, revised, and corrected the manuscript. All authors contributed to the article and approved the submitted version.

Funding

This work was provided by the Secretaría de Desarrollo Institucional Universidad Nacional Autónoma de México (UNAM) to EK-B and ID-L-C and from the CONACyT “Fronteras de la Ciencia” #1527 to JN-F. Postdoctoral research funds to EK-B were granted by the Direección General de Asuntos Academicos (DGAPA, UNAM).

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fevo.2022.916762/full#supplementary-material

Footnotes

1. ^ https://github.com/TransDecoder/TransDecoder
2. ^ https://solgenomics.net/
3. ^ http://evomics.org/wp-content/uploads/2016/06/cafe_tutorial.pdf

References