Chaoju Qian1
Xia Yan1,2
Hengxia Yin3
Xingke Fan1,4
Xiaoyue Yin1,4
Peipei Sun1,4
Zhijun Li5
Eviatar Nevo6*
Xiao-Fei Ma1*- 1Key Laboratory of Stress Physiology and Ecology in Cold and Arid Regions, Department of Ecology and Agriculture Research, Northwest Institute of Eco-Environment and Resources, Chinese Academy of Sciences, Lanzhou, China
- 2Key Laboratory of Ecohydrology of Inland River Basin, Northwest Institute of Eco-Environment and Resources, Chinese Academy of Sciences, Lanzhou, China
- 3State Key Laboratory of Plateau Ecology and Agriculture, Qinghai University, Xining, China
- 4College of Resources and Environment, University of Chinese Academy of Sciences, Beijing, China
- 5Xinjiang Production & Construction Corps Key Laboratory of Protection and Utilization of Biological Resources in Tarim Basin, Alar, China
- 6Institute of Evolution, University of Haifa, Haifa, Israel
As one of the hotspot regions for sympatric speciation studies, Evolution Canyon (EC) became an ideal place for its high level of microclimatic divergence interslopes. In this study, to highlight the genetic mechanisms of sympatric speciation, phenotypic variation on flowering time and transcriptomic divergence were investigated between two ecotypes of Ricotia lunaria, which inhabit the opposite temperate and tropical slopes of EC I (Lower Nahal Oren, Mount Carmel, Israel) separated by 100 m at the bottom of the slopes. Growth chamber results showed that flowering time of the ecotype from south-facing slope population # 3 (SFS 3) was significantly 3 months ahead of the north-facing slope population # 5 (NFS 5). At the same floral development stage, transcriptome analysis showed that 1,064 unigenes were differentially expressed between the two ecotypes, which enriched in the four main pathways involved in abiotic and/or biotic stresses responses, including flavonoid biosynthesis, α-linolenic acid metabolism, plant–pathogen interaction and linoleic acid metabolism. Furthermore, based on Ka/Ks analysis, nine genes were suggested to be involved in the ecological divergence between the two ecotypes, whose homologs functioned in RNA editing, ABA signaling, photoprotective response, chloroplasts protein-conducting channel, and carbohydrate metabolism in Arabidopsis thaliana. Among them, four genes, namely, SPDS1, FCLY, Tic21 and BGLU25, also showed adaptive divergence between R. lunaria and A. thaliana, suggesting that these genes could play an important role in plant speciation, at least in Brassicaceae. Based on results of both the phenotype of flowering time and comparative transcriptome, we hypothesize that, after long-time local adaptations to their interslope microclimatic environments, the molecular functions of these nine genes could have been diverged between the two ecotypes. They might differentially regulate the expression of the downstream genes and pathways that are involved in the interslope abiotic stresses, which could further diverge the flowering time between the two ecotypes, and finally induce the reproductive isolation establishment by natural selection overruling interslope gene flow, promoting sympatric speciation.
Introduction
Sympatric speciation, first proposed by Darwin as speciation occurring in contiguous populations with ongoing gene flow, was considered as one of the important models of biodiversity origin (Darwin, 1859). Compared to other geographic models of speciation, such as the allopatric and parapatric speciation models, sympatric speciation attracts particular interests for the reproductive barriers must evolve in situ to prevent homogenization, which was challenging to be proved between contiguous populations (Mayr, 1947, 1963; Fitzpatrick et al., 2008; Bird et al., 2012). Previous empirical studies and theoretical modeling supported that sympatric speciation was possible (Smith, 1966; Gavriletz, 2004; Bolnick and Fitzpatrick, 2007). However, it required the conditions of primary divergence with gene flow between populations, with complete absence of geographic barriers, which could be uncommon in nature. Recently, a plenty of studies on animal and bacteria revealed that the formation of reproductive isolation during the sympatric speciation could be favored by natural selection, such as in sexual selection and resource competence (Berlocher and Feder, 2002; Li K. et al., 2015, 2016; Rosser et al., 2015; Dunning et al., 2016; Getz et al., 2016; Kautt et al., 2016; Rodriguez et al., 2016; Zhao et al., 2016). Up to date, however, there is currently very little evidence to prove that natural selection in plants could induce sympatric speciation, partially due to the fact that flowering time hardly diverged between contiguous populations.
With high environmental heterogeneity, mountains supply diversifying habitats for biomes, which could have promoted the sympatric speciation. As one of the most promising models to address the sympatric speciation, “Evolution Canyon” (EC I) (Lower Nahal Oren, Mount Carmel, Israel) supplied plenty of cases in sympatric speciation across life from bacteria to mammals, which showed that similar genomic and phenotypic differentiation across species occurred between microsites with high microclimatical heterogeneity (Nevo, 2009). Such microclimatic differentiation has considerable consequences in population genetic structure within diverse organisms (Nevo, 1997, 2012; Rashkovetsky et al., 2006; Hadid et al., 2013, 2014), and the EC model reveals interslope divergence across life from bacteria to plants and mammals, followed by divergent genomic, proteomic, and phenomic adaptive complexes (Nevo, 2012, 2014, 2015; Hadid et al., 2014; Buse et al., 2015; Finkel et al., 2015). In the EC, major adaptive complexes on the tropical “African” slope [(AS), also called south-facing slope (SFS)] are related to solar radiation, heat, and drought, whereas those on the temperate “European” slope [(ES), also called north-facing slope (NFS)] are related to shade stress for photosynthesis, and interslope species divergence led to the incipient sympatric speciation in a diversity of organisms (Pavlícek et al., 2003; Nevo, 2006, 2014; Kossover et al., 2009; Yang et al., 2009; Li K. et al., 2016). This kind of local adaptation resulted from interslope microclimatic divergence also occurring in other three “Evolution Canyons” in Israel in the Galilee, Golan, and Negev Mountains as well as in the other evolution canyons worldwide (Nevo, 2012), and in the extension of the Evolution Canyon model in Evolution Plateau (upper Galilee Mountains) (Hadid et al., 2013; Li K. et al., 2015, 2016; Zhao et al., 2016).
Among the diverse organisms, an annual plant species of Brassicaceae (Warwick and Al-Shehbaz, 2006), Ricotia lunaria inhabits both contrasting slopes of EC I. Based on the AFLP (amplified fragment length polymorphism) and full genome array hybridizations in R. lunaria from the contrasting opposite slopes of EC I, previous studies showed that both the genetic diversity and transcriptome expression pattern differed significantly between interslope ecotypes. Up-regulation of drought resistance genes were witnessed in the SFS individuals, while the photosynthetic genes were upregulated in the NFS individuals, suggesting that the expression patterns diverged between interslopes due to environmental heterogeneity (Brodsky et al., 2008; Kossover et al., 2009). Unfortunately, little is known about how the microclimatic heterogeneity affected the phenotypic differentiation and even the genomic divergence between ecotypes from the two slopes. On the other hand, R. lunaria has relatively big and heavy seeds, and the gravity should play a significant role in its seed dispersal, suggesting higher intraslope than interslope gene flow. This supported the observation of higher genetic divergent interslopes than intraslopes (Kossover et al., 2009). For bottom slope populations at EC I, however, gravity did not intensify the gene flow between the slopes. Thus, the phenotypic and genetic divergence between the two bottom contrasting slope populations (#3 on SFS and #5 on NFS), could be better explained by microclimatic heterogeneity between the SFS and NFS.
In this study, we planted two ecotypes of R. lunaria from the bottom populations of the opposite slopes (population #3 on the dry tropical savannoid SFS and population #5 on the humid temperate forested NFS from EC I) (Figure 1). We investigated both the phenotypic variations of flowering time and the transcriptomic expression patterns between them. We aimed to address the following questions of adaptive interslopes micro-ecological divergence across ∼100 m: (1) whether the flowering times were significantly differentiated between the two interslope ecotypes; (2) whether the expression or sequence variations on the flowering time genes were also diverged between the two ecotypes? Likewise, what regulatory modules of genes or pathways that are adapted to the interslope microclimatic divergence induce the differentiation of flowering time between the two ecotypes? The investigation on both phenotypic and transcriptomic divergence of these two ecotypes of R. lunaria, will shed light on our better understanding of the molecular basis of sympatric speciation in plants in ecologically divergent microsites.
FIGURE 1. The opposing slopes of “Evolution Canyon” I, Lower Nahal Oren, Mount Carmel. The tropical xeric “African” slope is on the right, and the temperate mesic “European” slope is on the left. Panel (A) is the schematic diagram, panel (B) is the cross section view of EC I, Lower Nahal Oren, Mount Carmel, and panel (C) is the air view of EC I (Nevo, 2012), only one genotype from the station 5 and one genotype from the station 3 were analyzed in this study.
Materials and Methods
Plant Materials and Physiological Measurements
Two ecotypes of R. lunaria seeds were collected in the year 2012. As shown in Figure 1, one was from station 3 (SFS3) of EC I under savannoid tropical high illuminance, high temperature and dry growing conditions. By contrast, another one from the forested temperate opposite slope of EC I station 5 (NFS5), under the opposite growing conditions to the SFS3 ecotype with shade and humid environment. Before planting, the seeds were stored for 1 month at -20°C to synchronize their germinations. To verify whether genes under stress were constructively differently expressed between the two interslope ecotypes, five replicates for both the two ecotypes were planted in the same growth chamber at 23°C under long day conditions of 18 h light and 6 h dark, with 53% humidity conditions. The flowering times of the two ecotypes were recorded during the plants’ growth. Different tissues of the two ecotypes were collected for RNA isolation at the stage of the first flowering, and all materials were sampled at 10 am and immediately frozen in liquid nitrogen.
RNA Extraction, cDNA Library Preparation, and RNA-seq
Total RNA samples from root, stem, leaf, and bud for each ecotype were extracted using the E.Z.N.A. Plant RNA Kit (Omega Bio-tek, United States) with 2% PVBB added, and the remaining DNA was removed by RNase-free DNase (Omega Bio-tek, United States) according to the instruction manual. RNA concentration and purity of each sample was determined with ratio of OD260/280 by NanoDrop2000™ micro-volume spectrophotometer (Thermo Scientific, Waltham, MA, United States), and all the samples passed quality control as the ratio of OD260/280 between 1.9 and 2.2 and the ratio of OD260/230 was less than 2.0. Then, the RNA integrity was further verified by 1.5% agarose gel electrophoresis with two clear bands of 28S/18S ribosomal RNA. Total RNA samples for each ecotype was pooled together including the root, stem, leaf, and bud in the ratio of 1:1:1:2 to construct the cDNA library, respectively.
The poly-A mRNA was enriched with magnetic Oligo (dT) beads and then fragmented into short fragments as the templates to synthesize the cDNA fragments as detailed by Shi et al. (2013). The cDNA fragments were further purified with a QiaQuick PCR extraction kit (QIAGEN Inc., Valencia, CA, United States), followed by end repairing and tailing A, and then were ligated to sequencing adaptors. The required length ligation products with 200–700 bp were purified by agarose gel electrophoresis and enriched by PCR amplification. Finally, the paired-end library was sequenced by the Illumina HiSeqTM 2000 platform with the average read length of 90 bp.
De novo Assembly and Expression Profiling
After filtering adapter sequences, reads with ambigous sequences ‘N’ over 5% and reads with a base quality less than Q20, a total of 39.90 and 31.16 million clean reads were generated for NFS5 ecotype and SFS3 ecotype, respectively (Table 1). Then, Trinity was used to de novo assembly the high-quality clean reads into contigs (Grabherr et al., 2011). Short contigs were then clustered and assembled into the longest contigs with gaps according to the pair-end information and sequence similarity. The longest contigs in each cluster and singleton were combined together as the total unigenes. The RNA-seq raw data was deposited in the NCBI with SRA accession number of SRP150700.
To determine the unigene expression profiles, RPKM (reads per kilobases per million reads) of each unigene was computed by ERANGE3.1 software (Mortazavi et al., 2008). IDEG6 software (Romualdi et al., 2003) was also used to further identify the differentially expressed genes (DEGs) by pairwise comparison of the two ecotypes, with p-value < 0.01 and FDR significance score < 0.001 (Romualdi et al., 2003).
Gene Structure and Functional Annotation
The ORF (open reading frame) for each unigene was predicted by the Getorf program of the EMBOSS software package (Rice et al., 2000). SNPs between the two ecotypes were detected using SOAPsnp1 with the default parameters (Li et al., 2008). Functional annotation of each unigene was performed by aligning with public protein and/or nucleotide databases (such as the NCBI Nr, Nt databases, Swiss-Prot protein database, COG database, and the KEGG database) by BLASTx with an E-value cutoff of 1e-5. Further annotation analysis, such as the biological process, molecular functions and cellular components, were also performed with GO terms by Blast2GO software (Conesa et al., 2005; Conesa and Gotz, 2008).
Comparative Transcriptome Analysis and Putative Candidate Genes Involved in the Constructive Divergence of the Two Ecotypes
To further explore whether the constitute expression were divergent between the two ecotypes, differential expression analysis was first performed with EBSeq (Leng et al., 2013) under 1% false discover cutoff and at least twofold change. KEGG pathway analysis was also implemented to unravel the differential expression of genotype-dependent stress-responsive transcriptome pathways.
Based on SOAPsnp1, SNP calling was also performed between the transcriptomes of the two ecotypes to discover SNPs between the two genomes. Noteworthy, the whole-genome duplications in this genus (Lysak et al., 2007) could exaggerate the real number of the SNPs, as a large percentage of them could be due to the divergence between paralogous genes and unequal expression between the two ecotypes. Thus, the criteria were re-considered to guarantee as much as to further screen the SNPs on the orthologous genes, which were verified by the homologous unigene sequences based on very strict reciprocal localBlast, and then each of the unigene pairs were then trimmed into CDS fragments with SNPs.
We further identified the orthologous sequences with BLASTx against the unigenes library for each ecotype with Arabidopsis thaliana protein and/or nucleotide sequences from TAIR 11th public release2 under the threshold of E-value ≤ 1e-5. Furthermore, to screen whether the homologous genes diverged adaptively between the two ecotypes in response to their local environments between the slopes, intra-species (Pia/Pis) within the two ecotypes of R. lunaria and inter-species divergence (Ka/Ks) between R. lunaria and A. thaliana for each homologous genes were estimated with KaKs_Calculator 2.0 (Wang et al., 2010) under the YN model of approximate method (Yang and Nielsen, 2000).
Results and Discussion
Transcriptome Variations of the Two R. lunaria Ecotypes
Two transcriptomes with total RNA samples combined from bud, leaf, stem, and root (with the ratio of 2:1:1:1) were sequenced. In total, over 14 Gbp clean data from 71 million reads were sequenced and passed the Q20 quality control, with 8 Gbp for the ecotype of NFS5 and 6.2 Gbp for the ecotype of SFS3, respectively (Table 1). Interestingly, the SFS3 ecotype had slightly higher, 0.69%, GC content than that of the NFS5 ecotype (Table 1), which suggests that long-term higher solar radiation, heat, and drought stress could have induced the genome of R. lunaria to evolve adaptively to be more stable in the south-facing slope than the opposite slope. This was widely found in bacteria from the interslopes, such as Bacillus subtilis (Barash et al., 2006; Nevo, 2012), and Hallobacterium species (Kennedy et al., 2001).
TABLE 1. Information of transcriptomes and flowering time of the two genotypes from the contrasting slopes.
A total of 36,705 and 33,499 unigenes were assembled in the transcriptome of NFS5 ecotype and SFS3 ecotype, respectively. Both ecotype transcriptomes had N50 over 1,500 bp with the average length of 900 bp (Table 2). 47,103 unigenes were assembled when these two transcriptomes were pooled together, among which 14,341 (30.45%) unigenes were longer than 1 kb, suggesting the goodness of this de novo assembly. However, only 19,051 unigenes (40%) were shared with each other, while 16,163 unigenes were specific to NFS5 and 11,189 unigenes specific to SFS3. This phenomenon could have resulted for two reasons: (1) the one is for the vast divergence between the two ecotypes, which made it difficult to match the short reads together, since allelic SNPs and indels could introduce more isoforms by Trinity assembly (Grabherr et al., 2011); (2) the other one is for the incomplete sequencing of transcriptomes for both the two ecotypes.
TABLE 2. Length distribution of the two transcriptomes of R. lunaria.
Among the 47,103 unigenes, a total of 39,860 unigenes (84.6%, Table 3) could be annotated against with the public databases. Among these annotated unigenes, 34,878 unigenes (79.3%) can find the homologs in the Nt database with E-value ≤ 1e-5 (Table 3), and 10,559 unigenes out of them were classified into 24 COG categories (Figure 2A). Among these categories, the largest group is the “general function prediction only” (2,901, 27.47%), followed by “replication, recombination and repair” (1,511, 14.31%), “transcription” (1,449, 13.72%) and “signal transduction mechanisms” (1,266, 11.99%) (Figure 2A). A total of 32,191 unigenes were assigned into 10,942 GO terms of three main categories (cellular component, molecular function, and biological process) and 52 sub-categories (Figure 2B). Pathway-based analysis showed that 8,016 unigenes were assigned into 170 KEGG pathways, with the most highly represented category of “ribosome” (462 genes), “plant hormone signal transduction” (325 genes) and “oxidative phosphorylation” (287 genes) (Supplementary Table S1). Notably, most abundant genes in both of the two ecotypes were involved in oxidoreductase activity, photosynthesis, light harvesting, energy metabolism and defense responses to abiotic stress (Supplementary Table S2), which supported the hypothesis that environmental conditions in the growth chamber may be stressed for both ecotypes, especially for the light conditions, such as affecting population # 5 growing in shade.
FIGURE 2. Gene annotations of unigenes. (A) Histogram presentation of clusters of orthologous groups (COG) classification, (B) good hits were aligned to the GO database.
Constructive Differently Expressed Genes Between the Two Ecotypes in the Same Growth Chamber
To further explore whether the constitute expression was divergent between the two ecotypes, differential expression analysis was first performed with EBSeq (Leng et al., 2013) under 1% false discover cutoff and at least twofold change. Compared with the level of gene expression in SFS3 ecotype, a total of 1,064 unigenes were differentially expressed in the NFS5 ecotype, including 683 up-regulated and 381 down-regulated unigenes (Figure 3). Among the differentially expressed unigenes (DEGs), 553 unigenes were categorized in 21 COG clusters (E-value ≤ 1e-5) with the following five largest categories: general function prediction only (110 unigenes), carbohydrate transport and metabolism (51 unigenes), signal transduction mechanisms (49 unigenes), transcription (44 unigenes), and secondary metabolites biosynthesis, transport, and catabolism (38 unigenes) (Figure 4). Using the hypergeometric test compared to the genomic background, we further identified 42 significantly enriched GO terms of DEGs containing 976 unigenes (p ≤ 0.05, after Bonferroni correction, Supplementary Table S4), and all the top DEGs came from the pathways involved in biotic and abiotic responses, such as jasmonic acid biosynthetic process (GO:0009695), response to chitin (GO:0010200), response to fungus (GO:0009620), response to jasmonic acid stimulus (GO:0009753), response to mechanical stimulus (GO:0009612), abscisic acid mediated signaling pathway (GO:0009738), response to water deprivation (GO:0009414), and hyperosmotic salinity response (GO:0042538) (all the p-values as 0, Supplementary Table S4). Most of these processes were played a dominant role in stress defense for almost all the plant species analyzed, such as the jasmonic acid biosynthetic process (Robertseilaniantz et al., 2011; Pieterse et al., 2012).
FIGURE 3. Differentially expressed unigenes between NFS5 (T1) and SFS (T2). (A) Expression abundance comparison; (B) Volcano plot for the unigenes expressed differently between the two ecotypes. Red dots represented that unigenes with the expression level significantly higher in NFS5 (T1), while green dots for that of higher ones in SFS3 (T2).
FIGURE 4. Histogram presentation of clusters of orthologous groups (COG) classification for DEGs.
TABLE 3. Summary of gene annotation against the seven databases.
Differential Expression Pathways Analysis
Results of differential expression analysis implemented with KEGG pathway analysis showed that after genomic abundance correction (Supplementary Table S3), four pathways showed significantly different expression between the two ecotypes: flavonoid biosynthesis (10 unigenes, with 7 unigenes up-regulated in ecotype NFS5, while the other 3 unigenes up-regulated in ecotype SFS3, p = 7.13 E-06, Supplementary Figure S1), α-linolenic acid metabolism (14 unigenes, with 7 unigenes up-regulated in ecotype NFS5, and the other 7 unigenes up-regulated in ecotype SFS3, p = 1.77 E-05, Supplementary Figure S2), plant–pathogen interaction (26 unigenes, with 14 unigenes up-regulated in ecotype NFS5, while the other 12 unigenes up-regulated in ecotype SFS3, p = 2.59 E-04, Supplementary Figure S3) and linoleic acid metabolism (6 unigenes, with 3 unigenes up-regulated in ecotype NFS5, and another 3 unigenes up-regulated in ecotype SFS3, p = 1.3 E-02, Supplementary Figure S4). Flavonoid were reported as an important group of secondary metabolites, which played an important role in plant growth, development, and reproduction. It was also involved in the responses to diverse abiotic stresses, such as the UV radiation, high light, low availability of water and nutrients, temperature fluctuations and even the pathogen infection (Lenka et al., 2011; Ferreyra et al., 2012; Pandey et al., 2016; Singh et al., 2017). This kind of different expression of flavonoid biosynthesis pathway between different ecotypes were also found in the other model plants responding to diverse environment factors (Azuma et al., 2012; Guidi et al., 2016). In this study, several important enzymes in the flavonoid biosynthetic pathway were also differentially expressed between the two ecotypes (Supplementary Figure S1), suggesting that it could be an effective response to various stress-induced injuries in the growth chamber for R. lunaria, such as in light or shade stress for each ecotype. The two major unsaturated fatty acids in membrane lipids of plant leaves, α-linolenic and linoleic acid, also contribute to plant defense responses for abiotic and/or biotic stresses as documented in other plant species (Li Q. et al., 2016; Lim et al., 2017). The plant–pathogen interaction pathway was also reported to be not only responsive to the biotic stresses but also to the abiotic stresses in plants (Dodds and Rathjen, 2010; Barakate and Stephens, 2016; Wang et al., 2017). These four pathways were also widely involved in various plant resistance against diverse abiotic and/or biotic stress-induced injuries, such as in drought, UV, and/or salt responses for peach (Li X.W. et al., 2015), apple (Yin et al., 2016), rice (Lenka et al., 2011), Chorispora bungeana and Arabidopsis (Zhao et al., 2012). To respond to the ecological stresses in the growth chamber, R. lunaria could rely on the integrated effects of different pathways, and these up-regulated enzymes involved in these pathways could also provide targets for further studies of the molecular regulation mechanism in response to specific environmental stresses, such as the key enzymes of DFR, OPCL1, and JAZ. Clearly, the common garden experiment in the growth chamber would highlight ecotype divergence.
Genomic Divergence Between the Two Ecotypes
A total of 20,130 putative SNPs with scores over 30 were discovered in 6,107 genes. After filtering with strict reciprocal localBlast to avoid the paralogs, a total of 15,402 SNPs were finally found in 1,140 homologous genes between the two ecotypes, covering more than one million basepairs of the genome sequences (1.71 Mbp, with length varied from 164 bp to 6 kbp).
To further investigate whether the homologous genes diverged between the two ecotypes responding to the interslope divergent environments, among the 1,140 homologous genes with SNPs, we identified a set of 281 homologous genes with both synonymous and non-synonymous mutations (Figure 5). As shown in Table 4, nine pairs of homologous genes were putatively involved in the adaptive ecotypes’ divergence from each other with the Pia/Pis ratio over 1: T1_Unigene27281 vs. T2_Unigene22713 (RNA recognition motif-containing protein), T2_Unigene19717 vs. T1_Unigene31973 (WPP domain-interacting protein 3), T2_Unigene22075 vs. T1_Unigene22516 (spermidine synthase 1), T1_Unigene14300 vs. T2_Unigene15331 (regulatory component of ABA receptor 1), T1_Unigene30289 vs. T2_Unigene25344 (repressor of silencing 3), T1_Unigene25303 vs. T2_Unigene21895 (farnesylcysteine lyase), T1_Unigene24651 vs. T2_Unigene25310 (ZML2, GATA transcription factor 28), T1_Unigene24978 vs. T2_Unigene23451 (Translocon at inner membrane of chloroplasts 21), and T1_Unigene31300 vs. T2_Unigene16505 (BGLU25 beta glucosidase 25). Based on previous functional annotations, we found that most of these divergent genes were involved in the post-transcriptional regulation, ABA-dependent response, photoprotective response, chloroplast protein-conducting channel, carbohydrate metabolic, and so on. For examples, zinc finger protein expressed in inflorescence meristem like2 (ZML2), a GATA transcription factor, takes part in the photoprotective response mediated by the photoreceptor Cryptochrome 1 (cry1) (Shaikhali et al., 2012), which is a critical gene regulating flowering time in many plants (Yuan et al., 2016). ZML2 also plays a role in wound-induced lignification as a transcriptional repressor for the lignin biosynthesis in maize (Velezbermudez et al., 2015). The RNA recognition motif-containing (RRM), one of the largest families of RNA binding proteins gene widely detected not only in plants and animals, but also in fungus, plays an important role in RNA recognition and editing (Fetka et al., 2000; Rebay et al., 2000; Samarajeewa et al., 2017). In plants, as documented in model plants of Arabidopsis and maize, RRM is required for plastid and mitochondrial RNA editing (Sun et al., 2013), and previous studies showed that an organelle RRM mutant exhibits slower growth and delayed flowering time (Shi et al., 2016). As most of these divergent genes were located in the up-stream of the biological pathways, we speculated that the adaptive divergence could rely more on sensing of multiple stresses with upstream regulations than with downstream regulations. To verify this hypothesis, of course, more transgenic experiments on these divergent alleles are still needed.
FIGURE 5. Boxplot of the divergence on the 281 homologs within Ricotia lunaria (Pia, Pis, and Pia/Pis) and between R. lunaria and Arabidopsis thaliana (Ka, Ks, and Ka/Ks). ∗p-value < 0.05, ∗∗p-value < 0.01, which was estimated with Student’s t-test.
TABLE 4. Inter-species and intra-species divergence of genes involved into the ecological adaptation.
FIGURE 6. Simulated pathway of the ongoing sympatric speciation model in the two ecotypes of R. lunaria living on the opposite slopes of EC I.
Noteworthy, four of the isogenes (spermidine synthase 1, farnesylcysteine lyase, Translocon at inner membrane of chloroplasts 21, and BGLU25 beta glucosidase 25) also showed significant divergence between R. lunaria and A. thaliana with the Ka/Ks ratio over 1. spermidine synthase 1 (SPDS1), a member of the spermidine synthase-related gene family that takes part in the polyamine biosynthesis, is responsible for spermidine accumulation in many plants and expressed in all plant organs. It plays an important role in the regulation of spermidine synthase activity, and is involved in different growth stages by modulating the contents of polyamines in plant cells (Hanzawa et al., 2002; Neily et al., 2011). farnesylcysteine lyase (FCLY), encoding a specific farnesylcysteine (FC) lyase, expresses in all Arabidopsis tissues and organs. It is proved to be involved in the recycling and oxidative metabolism of FC to farnesal and cysteine, and in the regulation of ABA signaling and meristem development by regulating the accumulation of FC and inhibition of isoprenylcysteine methyltransferase (ICMT) (Crowell et al., 2007; Huizinga et al., 2010). translocon at inner membrane of chloroplasts 21 (Tic 21), also called permease in chloroplasts1 (PIC1), acts as an essential part in protein translocation across the inner envelope membrane of chloroplast and regulates plant growth and development by directing homeostasis and transports of iron. It also plays a critical role in leaf development in the later stages (Teng et al., 2006; Duy et al., 2007, 2011; Lopezmillan et al., 2016). beta glucosidase 25 (BGLU25) is involved in the carbohydrate metabolism and glycosyl compound metabolic process, which plays an important role in the regulation of hydrolyzing O-glycosyl compounds and activity of the beta-glucosidase and hydrolase (Xu et al., 2004). Considering that all of these four genes were also involved in the local adaptation between the two ecotypes of R. lunaria, we hypothesize that these genes could play important roles in plant speciation, at least in the Brassicaceae, and may be involved in the interslope sympatric speciation of R. lunaria in Evolution Canyon I. This phenomenon was also examined in Evolution Canyon II, as was done for Bacillus simplex (Sikorski and Nevo, 2005).
Phenotypic Variations and Regulation Differentiation of Flowering Time Between the Two Ecotypes of R. lunaria
Under long day conditions (with 18 h light with 6 h dark), we found out that the ecotype of SFS3, which is adapted to the dry tropical biome, is flowering at least 3 months earlier than that of the humid temperate ecotype NFS5, as evidenced by 57 days for ecotype SFS3, while 149 days for ecotype NFS5. This dramatic earlier phase transition from nutrition growth to reproductive growth under drought stress was widely found in many other natural plant species, such as in A. thaniana, Lupinus luteus L. and wheat (Franks et al., 2007; Berger and Ludwig, 2014; Farooq et al., 2014; Kenney et al., 2014). It was also documented that, to avoid the heat stress, plants tend to flower earlier (Suter and Widmer, 2013; Abouelwafa and Amein, 2016; Ishimaru et al., 2016). Nevo et al. (2012) showed that ∼10 days earlier flowering in 10 wild barley, Hordeum spontaneum, and 10 wild emmer wheat, Triticum dicoccoides, populations across Israel, due to global warming between 1980 and 2008, highly support that earlier flowering is an effective way for plants to escape from both the drought and heat stresses.
Previous studies displayed by tiling array hybridizations on the wild sampling from the opposite slopes of EC I showed that DEGs were involved in at least the blue light signaling pathway, circadian rhythm and protein amino acid phosphorylation in response to the drought stress on the SFS. Interestingly, more heat-response genes are activated and a fivefold increase in flower development was recorded in the SFS ecotypes of R. lunaria (Brodsky et al., 2008). Although the circadian rhythm – plant pathway is not significantly differently expressed between the two ecotypes after genomic abundance correction, four up-regulated genes encoding three key enzymes were still detected between them, namely, T1_Unigene26319 (APR9); T1_Unigene31423 and T2_Unigene30678 (LHY); T2_Unigene24716 (CCA1) (Supplementary Figure S5). Until now, functions of these core circadian clock genes in plant’ stress response are still unclear. In this study, our controlled common garden experiments in the growth chamber witnessed both huge phenotypic (mainly on flowering time) and genomic divergence between the two interslope ecotypes. From our data, we found out that nine adaptive genes were involved in the post-transcriptional regulation, ABA-dependent response, photoprotective response, chloroplast protein-conducting channel and carbohydrate metabolic pathways, which might play a role in sensing and responding to the aforementioned multiple stresses. The huge divergence on their sequences between the two ecotypes indicated that they might have been also diversified with allelic functions, which were favored by the microecologically climatic heterogeneity of the interslopes during the sympatric speciation.
Considering that the expression patterns were different between the two ecotypes, both in wild and growth chamber conditions, the two ecotypes might have evolved divergent regulation systems responding to multiple stresses. In the model plant Arabidopsis, the flowering time was orchestrated by many pathways, including the photoperiodic pathway, vernalization pathway, autonomous pathway, and gibberellins (GA) pathway (Boss et al., 2004; Mulekar and Huq, 2012; Osnato et al., 2012; Kemi et al., 2013). As a close relative species of Arabidopsis, the regulation of flowering time in R. lunaria could be similar to the Arabidopsis species. Although plenty of studies supported that drought and heat stresses could promote the flowering time in plants (Franks et al., 2007; Suter and Widmer, 2013; Berger and Ludwig, 2014; Farooq et al., 2014; Kenney et al., 2014; Abouelwafa and Amein, 2016; Ishimaru et al., 2016), the crosstalk between modules is still unclear. In this study, we found out that the ZML2 gene was under divergence selection, which could be integrated into the photoreceptor CRY1 responding to the photoprotection (Shaikhali et al., 2012). As CRY1 is a well-known upstream geen involved in flowering time regulation (Yuan et al., 2016), we propose that ZML2 might also play an important role in the differentiation of the flowering time between the two ecotypes as a consequence of local interslope adaptation. Thus, there might be some new pathways involved in the regulation of flowering time in response to the multiple stresses (Figure 6). After long-term local adaptation to the interslope microclimatically contrasting environments, both the phenotype traits and genomes of these two ecotypes of R. lunaria have already greatly diverged from each other, and further induced the differentiation of expression patterns of numerous genes. Under the same common garden environmental conditions in the growth chamber, the stresses would be different for each ecotype. For example, the conditions in the growth chamber could be shade stress for the SFS3 ecotype, while it could be high light stress for the NFS5 ecotype. To respond to multiple stresses, the common genes and pathways involved in biotic and abiotic stresses significantly differentially expressed multiple genes between the two ecotypes, which might consequently impact the expression patterns of the down-stream flowering time genes such as FT1 (FLOWERING LOCUS T1). As a long evolutionary outcome, the flowering time of the two ecotypes have dramatically diverged, and further developed reproductive isolation. Due to the diversifying selection overruling gene flow (Nevo, 2011), this kind of incipient sympatrical speciation between the two ecotypes in Evolution Canyon I was also found in other species, such as wild barley (Nevo, 2006) and Drosophila melanogaster (Korol et al., 2006; Pavlíèek et al., 2008; Nevo, 2011). Of course, more genetic evidences are still needed to verify the molecular functions of these adaptive alleles in R. lunaria.
Conclusion
Sympatric speciation is a common model widely existing across the canyons with climatic, geologic, and biotic differentiation between slopes. In this study, we fully investigated a new model for the incipient sympatric speciation at Evolution Canyon (Nevo, 2011). Based on the phenotypic investigation of flowering time in the common garden, we identified 3 months earlier flowering time of the dry tropical sun population of R. lunaria than the humid-cool temperate shade population. Likewise, between the opposite slope ecotypes, a total of 1,064 isogenic genes were significantly differentially expressed and nine genes significantly diverged that displayed specific interslope genetic-ecological adaptations. Both the flowering time and the expression patterns of the genes related to multiple stresses significantly differed between the two ecotypes. Clearly, flowering time is a cardinal variable of reproductive isolation, as a consequence of the genomic divergence between plant populations. This work substantiated a new species to the increasing list of either incipiently or fully sympatrically evolving species in Evolution Canyon.
As the other models in Evolution Canyon (Sikorski and Nevo, 2005; Nevo, 1995, 2014). Evolution Plateau (Hadid et al., 2013; Li K. et al., 2015, 2016; Zhao et al., 2016) and Evolution Slope (Wang et al., 2018), sympatric speciation of R. lunaria at Evolution Canyon I in Mount Carmel could also be tested at Evolution II in western Upper Galilee, which are separated by 40 km from EC I. Future projects of R. lunaria at both EC I and EC II could include intraslope and interslope crosses, investigations on genomics, transcriptomics, epigenomics, and metabolomics, which could shed light on the complexities of sympatric speciation.
Author Contributions
X-FM and EN conceived the ideas and designed the investigations. CQ, XYa, HY, and X-FM analyzed the data and wrote the first version of the manuscript. CQ and HY collected the data. XF, XYa, and PS assisted in performing the experiments. ZL and EN provided improvements to the manuscript.
Funding
This work was supported by grants from the National Natural Science Foundation of China (NSFC, Grant Nos. 31370395 and 31500266), the “One Hundred Talents” project of the Chinese Academy of Sciences (Grant No. 29Y127E71), the “Ancell-Teicher Research Foundation of Genetics and Molecular Evolution” of Israel, and the Xinjiang Production and Construction Crops Key Laboratory of Protection and Utilization of Biological Resources in Tarim Basin (Grant No. BRZD1809).
Conflict of Interest Statement
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.
Acknowledgments
We thank Xin Zhao and Pengshan Zhao for laboratory assistance.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2018.00506/full#supplementary-material
Footnotes
References
Abouelwafa, S. F., and Amein, K. A. (2016). Genetic diversity and potential high temperature tolerance in barley (Hordeum vulgare). World J. Agric. Res. 4, 1–8.
Azuma, A., Yakushiji, H., Koshita, Y., and Kobayashi, S. (2012). Flavonoid biosynthesis-related genes in grape skin are differentially regulated by temperature and light conditions. Planta 236, 1067–1080. doi: 10.1007/s00425-012-1650-x
Barakate, A., and Stephens, J. (2016). An overview of CRISPR-based tools and their improvements: new opportunities in understanding plant–pathogen interactions for better crop protection. Front. Plant Sci. 7:765. doi: 10.3389/fpls.2016.00765
Barash, D., Sikorski, J., Perry, E. B., Nevo, E., and Nudler, E. (2006). Adaptive mutations in RNA-based regulatory mechanisms: computational and experimental investigations. Isr. J. Ecol. Evol. 52, 263–279. doi: 10.1560/IJEE_52_3-4_263
Berger, J., and Ludwig, C. (2014). Contrasting adaptive strategies to terminal drought-stress gradients in Mediterranean legumes: phenology, productivity, and water relations in wild and domesticated Lupinus luteus L. J. Exp. Bot. 65, 6219–6229. doi: 10.1093/jxb/eru006
Berlocher, S. H., and Feder, J. L. (2002). Sympatric speciation in phytophagous insects: moving beyond controversy? Annu. Rev. Entomol. 47, 773–815. doi: 10.1146/annurev.ento.47.091201.145312
Bird, C. E., Fernandez-Silva, I., Skillings, D. J., and Toonen, R. J. (2012). Sympatric speciation in the post “modern synthesis” era of evolutionary biology. Evol. Biol. 39, 158–180. doi: 10.1007/s11692-012-9183-6
Bolnick, D. I., and Fitzpatrick, B. M. (2007). Sympatric speciation: models and empirical evidence. Annu. Rev. Ecol. Evol. Syst. 38, 459–487. doi: 10.1146/annurev.ecolsys.38.091206.095804
Boss, P. K., Bastow, R., Mylne, J. S., and Dean, C. (2004). Multiple pathways in the decision to flower: enabling, promoting, and resetting. Plant Cell 16(Suppl.), S18–S31. doi: 10.1105/tpc.015958
Brodsky, L., Kossover, O., Ben-Jacob, E., and Korol, A. (2008). “Genome transcription divergent patterns of Ricotia lunaria on opposite slopes of “Evolution Canyon” revealed by tiling array hybridizations,” in Proceedings of the Plant and Animal Genomes XVI Conference, San Diego, CA, 12–16.
Buse, J., Fassbender, S., Entling, M. H., and Pavlicek, T. (2015). Microclimatic divergence in a mediterranean canyon affects richness, composition, and body size in saproxylic beetle assemblages. PLoS One 10:e0129323. doi: 10.1371/journal.pone.0129323
Conesa, A., and Gotz, S. (2008). Blast2GO: a comprehensive suite for functional analysis in plant genomics. Int. J. Plant Genomics 2008:619832. doi: 10.1155/2008/619832
Conesa, A., Gotz, S., Garcia-Gomez, J. M., Terol, J., Talon, M., and Robles, M. (2005). Blast2GO: a universal tool for annotation, visualization and analysis in functional genomics research. Bioinformatics 21, 3674–3676. doi: 10.1093/bioinformatics/bti610
Crowell, D. N., Huizinga, D. H., Deem, A., Trobaugh, C. W., Denton, R., and Sen, S. E. (2007). Arabidopsis thaliana plants possess a specific farnesylcysteine lyase that is involved in detoxification and recycling of farnesylcysteine. Plant J. 50, 839–847. doi: 10.1111/j.1365-313X.2007.03091.x
Dodds, P. N., and Rathjen, J. P. (2010). Plant immunity: towards an integrated view of plant–pathogen interactions. Nat. Rev. Genet. 11, 539–548. doi: 10.1038/nrg2812
Dunning, L. T., Hipperson, H., Baker, W. J., Butlin, R. K., Devaux, C., Hutton, I., et al. (2016). Ecological speciation in sympatric palms: 1. Gene expression, selection and pleiotropy. J. Evol. Biol. 29, 1472–1487. doi: 10.1111/jeb.12895
Duy, D., Stube, R., Wanner, G., and Philippar, K. (2011). The chloroplast permease PIC1 regulates plant growth and development by directing homeostasis and transport of iron. Plant Physiol. 155, 1709–1722. doi: 10.1104/pp.110.170233
Duy, D., Wanner, G., Meda, A. R., Von Wiren, N., Soll, J., and Philippar, K. (2007). PIC1, an ancient permease in Arabidopsis chloroplasts, mediates iron transport. Plant Cell 19, 986–1006. doi: 10.1105/tpc.106.047407
Farooq, M., Hussain, M., and Siddique, K. H. (2014). Drought stress in wheat during flowering and grain-filling periods. Crit. Rev. Plant Sci. 33, 331–349. doi: 10.3389/fpls.2017.01950
Ferreyra, M. L. F., Rius, S. P., and Casati, P. (2012). Flavonoids: biosynthesis, biological functions, and biotechnological applications. Front. Plant Sci. 3:222. doi: 10.3389/fpls.2012.00222
Fetka, I., Radeghieri, A., and Bouwmeester, T. (2000). Expression of the RNA recognition motif-containing protein SEB-4 during Xenopus embryonic development. Mech. Dev. 94, 283–286. doi: 10.1016/S0925-4773(00)00284-7
Finkel, M., Ofer, J., Beharav, A., and Nevo, E. (2015). Species interslope divergence of ants caused by sharp microclimatic stresses at ‘Evolution Canyon’ II, Lower Nahal Keziv, western Upper Galilee, Israel. Isr. J. Entomol. 44, 63–73.
Fitzpatrick, B. M., Fordyce, J. A., and Gavrilets, S. (2008). What, if anything, is sympatric speciation? J. Evol. Biol. 21, 1452–1459. doi: 10.1111/j.1420-9101.2008.01611.x
Franks, S. J., Sim, S., and Weis, A. E. (2007). Rapid evolution of flowering time by an annual plant in response to a climate fluctuation. Proc. Natl. Acad. Sci. U.S.A. 104, 1278–1282. doi: 10.1073/pnas.0608379104
Gavriletz, S. (2004). Fitness Landscapes and the Origin of Species (MPB-41). Princeton, NJ: Princeton University Press.
Getz, W. M., Salter, R. M., Seidel, D. P., and Van Hooft, P. (2016). Sympatric speciation in structureless environments. BMC Evol. Biol. 16:50. doi: 10.1186/s12862-016-0617-0
Grabherr, M. G., Haas, B. J., Yassour, M., Levin, J. Z., Thompson, D. A., Amit, I., et al. (2011). Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nat. Biotechnol. 29, 644–652. doi: 10.1038/nbt.1883
Guidi, L., Brunetti, C., Fini, A., Agati, G., Ferrini, F., Gori, A., et al. (2016). UV radiation promotes flavonoid biosynthesis, while negatively affecting the biosynthesis and the de-epoxidation of xanthophylls: consequence for photoprotection? Environ. Exp. Bot. 127, 14–25. doi: 10.1016/j.envexpbot.2016.03.002
Hadid, Y., Pavlíèek, T., Beiles, A., Ianovici, R., Raz, S., and Nevo, E. (2014). Sympatric incipient speciation of spiny mice Acomys at “Evolution Canyon,” Israel. Proc. Natl. Acad. Sci. U.S.A. 111, 1043–1048. doi: 10.1073/pnas.1322301111
Hadid, Y., Tzur, S., Pavlíèek, T., Šumbera, R., Šklíba, J., Lövy, M., et al. (2013). Possible incipient sympatric ecological speciation in blind mole rats (Spalax). Proc. Natl. Acad. Sci. U.S.A. 110, 2587–2592. doi: 10.1073/pnas.1222588110
Hanzawa, Y., Imai, A., Michael, A. J., Komeda, Y., and Takahashi, T. (2002). Characterization of the spermidine synthase-related gene family in Arabidopsis thaliana. FEBS Lett. 527, 176–180. doi: 10.1016/S0014-5793(02)03217-9
Huizinga, D. H., Denton, R., Koehler, K. G., Tomasello, A., Wood, L., Sen, S. E., et al. (2010). Farnesylcysteine lyase is involved in negative regulation of abscisic acid signaling in Arabidopsis. Mol. Plant 3, 143–155. doi: 10.1093/mp/ssp091
Ishimaru, T., Hirabayashi, H., Sasaki, K., Ye, C., and Kobayashi, A. (2016). Breeding efforts to mitigate damage by heat stress to spikelet sterility and grain quality. Plant Prod. Sci. 19, 12–21. doi: 10.1080/1343943X.2015.1128113
Kautt, A. F., Machadoschiaffino, G., and Meyer, A. (2016). Multispecies outcomes of sympatric speciation after admixture with the source population in two radiations of Nicaraguan crater lake cichlids. PLoS Genet. 12:e1006157. doi: 10.1371/journal.pgen.1006157
Kemi, U., Niittyvuopio, A., Toivainen, T., Pasanen, A., Quilotturion, B., Holm, K., et al. (2013). Role of vernalization and of duplicated FLOWERING LOCUS C in the perennial Arabidopsis lyrata. New Phytol. 197, 323–335. doi: 10.1111/j.1469-8137.2012.04378.x
Kennedy, S. P., Ng, W. V., Salzberg, S. L., Hood, L., and DasSarma, S. (2001). Understanding the adaptation of Halobacterium species NRC-1 to its extreme environment through computational analysis of its genome sequence. Genome Res. 11, 1641–1650. doi: 10.1101/gr.190201
Kenney, A. M., McKay, J. K., Richards, J. H., and Juenger, T. E. (2014). Direct and indirect selection on flowering time, water-use efficiency (WUE, δ13C), and WUE plasticity to drought in Arabidopsis thaliana. Ecol. Evol. 4, 4505–4521. doi: 10.1002/ece3.1270
Korol, A., Rashkovetsky, E., Iliadi, K., and Nevo, E. (2006). Drosophila flies in “Evolution Canyon” as a model for incipient sympatric speciation. Proc. Natl. Acad. Sci. U.S.A. 103, 18184–18189. doi: 10.1073/pnas.0608777103
Kossover, O., Frenkel, Z., Korol, A., and Nevo, E. (2009). Genetic diversity and stress of Ricotia lunaria in “Evolution Canyon,” Israel. J. Hered. 100, 432–440. doi: 10.1093/jhered/esp014
Leng, N., Dawson, J. A., Thomson, J. A., Ruotti, V., Rissman, A. I., Smits, B. M., et al. (2013). EBSeq: an empirical Bayes hierarchical model for inference in RNA-seq experiments. Bioinformatics 29, 1035–1043. doi: 10.1093/bioinformatics/btt087
Lenka, S. K., Katiyar, A., Chinnusamy, V., and Bansal, K. C. (2011). Comparative analysis of drought-responsive transcriptome in Indica rice genotypes with contrasting drought tolerance. Plant Biotechnol. J. 9, 315–327. doi: 10.1111/j.1467-7652.2010.00560.x
Li, K., Hong, W., Jiao, H., Wang, G. D., Rodriguez, K. A., Buffenstein, R., et al. (2015). Sympatric speciation revealed by genome-wide divergence in the blind mole rat Spalax. Proc. Natl. Acad. Sci. U.S.A. 112, 11905–11910. doi: 10.1073/pnas.1514896112
Li, X. W., Jiang, J., Zhang, L., Yu, Y., Ye, Z., Wang, X., et al. (2015). Identification of volatile and softening-related genes using digital gene expression profiles in melting peach. Tree Genet. Genomes 11:71. doi: 10.1007/s11295-015-0891-9
Li, K., Wang, H., Cai, Z., Wang, L., Xu, Q., Lovy, M., et al. (2016). Sympatric speciation of spiny mice, Acomys, unfolded transcriptomically at Evolution Canyon, Israel. Proc. Natl. Acad. Sci. U.S.A. 113, 8254–8259. doi: 10.1073/pnas.1608743113
Li, Q., Lei, S., Du, K., Li, L., Pang, X., Wang, Z., et al. (2016). RNA-seq based transcriptomic analysis uncovers α-linolenic acid and jasmonic acid biosynthesis pathways respond to cold acclimation in Camellia japonica. Sci. Rep. 6:36463. doi: 10.1038/srep36463
Li, R. Q., Li, Y. R., Kristiansen, K., and Wang, J. (2008). SOAP: short oligonucleotide alignment program. Bioinformatics 24, 713–714. doi: 10.1093/bioinformatics/btn025
Lim, G., Singhal, R., Kachroo, A., and Kachroo, P. (2017). Fatty acid– and lipid-mediated signaling in plant defense. Annu. Rev. Phytopathol. 55, 505–536. doi: 10.1146/annurev-phyto-080516-035406
Lopezmillan, A. F., Duy, D., and Philippar, K. (2016). Chloroplast iron transport proteins – function and impact on plant physiology. Front. Plant Sci. 7:178. doi: 10.3389/fpls.2016.00178
Lysak, M. A., Cheung, K., Kitschke, M., and Bureš, P. (2007). Ancestral chromosomal blocks are triplicated in Brassiceae species with varying chromosome number and genome size. Plant Physiol. 145, 402–410. doi: 10.1104/pp.107.104380
Mayr, E. (1947). Systematics and the Origin of Species. Washington, DC: The National Academies Press.
Mayr, E. (1963). Animal Speciation and Evolution. Cambridge, MA: Har-vard University Press. doi: 10.4159/harvard.9780674865327
Mortazavi, A., Williams, B. A., Mccue, K., Schaeffer, L., and Wold, B. (2008). Mapping and quantifying mammalian transcriptomes by RNA-Seq. Nat. Methods 5, 621–628. doi: 10.1038/nmeth.1226
Mulekar, J. J., and Huq, E. (2012). Does CK2 affect flowering time by modulating the autonomous pathway in Arabidopsis. Plant Signal. Behav. 7, 292–294. doi: 10.4161/psb.18883
Neily, M. H., Baldet, P., Arfaoui, I., Saito, T., Li, Q., Asamizu, E., et al. (2011). Overexpression of apple spermidine synthase 1 (MdSPDS1) leads to significant salt tolerance in tomato plants. Plant Biotechnol. 28, 33–42. doi: 10.5511/plantbiotechnology.10.1013a
Nevo, E. (1995). Asian, African and European biota meet at ‘Evolution Canyon’ Israel: local tests of global biodiversity and genetic diversity patterns. Proc. R. Soc. B Bio. Sci. 262, 149–155. doi: 10.1098/rspb.1995.0189
Nevo, E. (1997). Evolution in action across phylogeny caused by microclimatic stresses at “Evolution Canyon”. Theor. Popul. Biol. 52, 231–243. doi: 10.1006/tpbi.1997.1330
Nevo, E. (2006). ” Evolution Canyon”: a microcosm of life’s evolution focusing on adaptation and speciation. Isr. J. Ecol. Evol. 52, 501–506. doi: 10.1560/IJEE_52_3-4_485
Nevo, E. (2009). Evolution in action across life at “Evolution Canyons”, Israel. Trends Evol. Biol. 1:3. doi: 10.4081/eb.2009.e3
Nevo, E. (2011). “Selection overrules gene flow at “Evolution Canyons”, Israel,” in Advance in Genetics Research, Vol. 5, ed. K. V. Urbano (New York, NY: Nova Science Publishers Inc.) 67–89.
Nevo, E. (2012). “Evolution Canyon,” a potential microscale monitor of global warming across life. Proc. Natl. Acad. Sci. U.S.A. 109, 2960–2965. doi: 10.1073/pnas.1120633109
Nevo, E. (2014). Evolution in action: adaptation and incipient sympatric speciation with gene flow across life at “Evolution Canyon”, Israel. Isr. J. Ecol. Evol. 60, 85–98. doi: 10.1080/15659801.2014.986879
Nevo, E. (2015). Evolution of wild barley at “Evolution Canyon”: adaptation, speciation, pre-agricultural collection, and barley improvement. Isr. J. Plant Sci. 62, 22–32. doi: 10.1080/07929978.2014.940783
Nevo, E., Fu, Y. B., Pavlicek, T., Khalifa, S., Tavasi, M., and Beiles, A. (2012). Evolution of wild cereals during 28 years of global warming in Israel. Proc. Natl. Acad. Sci. U.S.A. 109, 3412–3415. doi: 10.1073/pnas.1121411109
Osnato, M., Castillejo, C., Matiashernandez, L., and Pelaz, S. (2012). TEMPRANILLO genes link photoperiod and gibberellin pathways to control flowering in Arabidopsis. Nat. Commun. 3:808. doi: 10.1038/ncomms1810
Pandey, A., Alok, A., Lakhwani, D., Singh, J., Asif, M. H., and Trivedi, P. K. (2016). Genome-wide expression analysis and metabolite profiling elucidate transcriptional regulation of flavonoid biosynthesis and modulation under abiotic stresses in banana. Sci. Rep. 6:31361. doi: 10.1038/srep31361
Pavlíèek, T., Frenkel, Z., Korol, A. B., Beiles, A., and Nevo, E. (2008). Drosophila at the” Evolution Canyon” microsite, Mt. Carmel, Israel: selection overrules migration. Isr. J. Ecol. Evol. 54, 165–180. doi: 10.1080/15659801.2008.10639620
Pavlícek, T., Sharon, D., Kravchenko, V., Saaroni, H., and Nevo, E. (2003). Microclimatic interslope differences underlying biodiversity contrasts in” Evolution Canyon”, Mt. Carmel, Israel. Isr. J. Earth Sci. 52, 1–9. doi: 10.1560/QD0A-9T22-4FGL-FGH6
Pieterse, C. M. J., Der Does, D. V., Zamioudis, C., Leonreyes, A., and Van Wees, S. C. M. (2012). Hormonal modulation of plant immunity. Annu. Rev. Cell Dev. Biol. 28, 489–521. doi: 10.1146/annurev-cellbio-092910-154055
Rashkovetsky, E., Iliadi, K., Michalak, P., Lupu, A., Nevo, E., Feder, M., et al. (2006). Adaptive differentiation of thermotolerance in Drosophila along a microclimatic gradient. Heredity 96, 353–359. doi: 10.1038/sj.hdy.6800784
Rebay, I., Chen, F., Hsiao, F., Kolodziej, P. A., Kuang, B. H., Laverty, T. R., et al. (2000). A genetic screen for novel components of the Ras/Mitogen-activated protein kinase signaling pathway that interact with the yan gene of Drosophila identifies split ends, a new RNA recognition motif-containing protein. Genetics 154, 695–712.
Rice, P., Longden, I., and Bleasby, A. (2000). EMBOSS: the European molecular biology open software suite. Trends Genet. 16, 276–277. doi: 10.1016/S0168-9525(00)02024-2
Robertseilaniantz, A., Grant, M., and Jones, J. D. G. (2011). Hormone crosstalk in plant disease and defense : more than just JASMONATE-SALICYLATE antagonism. Annu. Rev. Phytopathol. 49, 317–343. doi: 10.1146/annurev-phyto-073009-114447
Rodriguez, K. A., Li, K., Nevo, E., and Buffenstein, R. (2016). Mechanisms regulating proteostasis are involved in sympatric speciation of the blind mole rat, Spalax galili. Autophagy 12, 703–704. doi: 10.1080/15548627.2016.1143592
Romualdi, C., Bortoluzzi, S., D’Alessi, F., and Danieli, G. A. (2003). IDEG6: a web tool for detection of differentially expressed genes in multiple tag sampling experiments. Physiol. Genomics 12, 159–162. doi: 10.1152/physiolgenomics.00096.2002
Rosser, N., Kozak, K. M., Phillimore, A. B., and Mallet, J. (2015). Extensive range overlap between heliconiine sister species: evidence for sympatric speciation in butterflies? BMC Evol. Biol. 15:125. doi: 10.1186/s12862-015-0420-3
Samarajeewa, D. A., Manitchotpisit, P., Henderson, M., Xiao, H., Rehard, D. G., Edwards, K. A., et al. (2017). An RNA recognition motif-containing protein functions in meiotic silencing by unpaired DNA. G3 7, 2871–2882. doi: 10.1534/g3.117.041848
Shaikhali, J., Barajaslopez, J. D. D., Otvos, K., Kremnev, D., Garcia, A. S., Srivastava, V., et al. (2012). The CRYPTOCHROME1-dependent response to excess light is mediated through the transcriptional activators ZINC FINGER PROTEIN EXPRESSED IN INFLORESCENCE MERISTEM LIKE1 and ZML2 in Arabidopsis. Plant Cell 24, 3009–3025. doi: 10.1105/tpc.112.100099
Shi, X., Germain, A., Hanson, M. R., and Bentolila, S. (2016). RNA recognition motif-containing protein ORRM4 broadly affects mitochondrial RNA editing and impacts plant development and flowering. Plant Physiol. 170, 294–309. doi: 10.1104/pp.15.01280
Shi, Y., Yan, X., Zhao, P. S., Yin, H. X., Zhao, X., Xiao, H. L., et al. (2013). Transcriptomic analysis of a tertiary relict plant, extreme xerophyte Reaumuria soongorica to identify genes related to drought adaptation. PLoS One 8:e63993. doi: 10.1371/journal.pone.0063993
Sikorski, J., and Nevo, E. (2005). Adaptation and incipient sympatric speciation of Bacillus simplex under microclimatic contrast at “Evolution Canyons” I and II, Israel. Proc. Natl. Acad. Sci. U.S.A. 102, 15924–15929. doi: 10.1073/pnas.0507944102
Singh, O., Kaur, R., and Mahajan, R. K. (2017). Flavonoid-surfactant interactions: a detailed physicochemical study. Spectrochim. Acta A Mol. Biomol. Spectrosc. 170, 77–88. doi: 10.1016/j.saa.2016.07.007
Sun, T., Germain, A., Giloteaux, L., Hammani, K., Barkan, A., Hanson, M. R., et al. (2013). An RNA recognition motif-containing protein is required for plastid RNA editing in Arabidopsis and maize. Proc. Natl. Acad. Sci. U.S.A. 110, E1169–E1178. doi: 10.1073/pnas.1220162110
Suter, L., and Widmer, A. (2013). Phenotypic effects of salt and heat stress over three generations in Arabidopsis thaliana. PLoS One 8:e80819. doi: 10.1371/journal.pone.0080819
Teng, Y., Su, Y., Chen, L., Lee, Y., Hwang, I., and Li, H. (2006). Tic21 is an essential translocon component for protein translocation across the chloroplast inner envelope membrane. Plant Cell 18, 2247–2257. doi: 10.1105/tpc.106.044305
Velezbermudez, I., Salazarhenao, J. E., Fornale, S., Lopezvidriero, I., Francozorrilla, J., Grotewold, E., et al. (2015). A MYB/ZML complex regulates wound-induced lignin genes in maize. Plant Cell 27, 3245–3259. doi: 10.1105/tpc.15.00545
Wang, D., Zhang, Y., Zhang, Z., Zhu, J., and Yu, J. (2010). KaKs_Calculator 2.0: a toolkit incorporating gamma-series methods and sliding window strategies. Genomics Proteomics Bioinformatics 8, 77–80. doi: 10.1016/S1672-0229(10)60008-3
Wang, M., Gao, L., Dong, S., Sun, Y., Shen, Q., and Guo, S. (2017). Role of silicon on plant–pathogen interactions. Front. Plant Sci. 8:701. doi: 10.3389/fpls.2017.00701
Wang, X., Chen, Z. H., Yang, C., Zhang, X., Jin, G., Chen, G., et al. (2018). Genomic adaptation to drought in wild barley is driven by edaphic natural selection at the Tabigha Evolution Slope. Proc. Natl. Acad. Sci. U.S.A. 115, 5223–5228. doi: 10.1073/pnas.1721749115
Warwick, S., and Al-Shehbaz, I. (2006). Brassicaceae: chromosome number index and database on CD-Rom. Plant Syst. Evol. 259, 237–248. doi: 10.1007/s00606-006-0421-1
Xu, Z., Escamillatrevino, L. L., Zeng, L., Lalgondar, M., Bevan, D. R., Winkel, B. S. J., et al. (2004). Functional genomic analysis of Arabidopsis thaliana glycoside hydrolase family 1. Plant Mol. Biol. 55, 343–367. doi: 10.1007/s11103-004-0790-1
Yang, Z., and Nielsen, R. (2000). Estimating synonymous and nonsynonymous substitution rates under realistic evolutionary models. Mol. Biol. Evol. 17, 32–43. doi: 10.1093/oxfordjournals.molbev.a026236
Yang, Z., Zhang, T., Bolshoy, A., Beharav, A., and Nevo, E. (2009). Adaptive microclimatic structural and expressional dehydrin 1 evolution in wild barley, Hordeum spontaneum, at ‘Evolution Canyon’, Mount Carmel, Israel. Mol. Ecol. 18, 2063–2075. doi: 10.1111/j.1365-294X.2009.04140.x
Yin, Z., Ke, X., Kang, Z., and Huang, L. (2016). Apple resistance responses against Valsa mali revealed by transcriptomics analyses. Physiol. Mol. Plant Pathol. 93, 85–92. doi: 10.1016/j.pmpp.2016.01.004
Yuan, S., Zhang, Z., Zheng, C., Zhao, Z., Wang, Y., Feng, L., et al. (2016). Arabidopsis cryptochrome 1 functions in nitrogen regulation of flowering. Proc. Natl. Acad. Sci. U.S.A. 113, 7661–7666. doi: 10.1073/pnas.1602004113
Zhao, Y., Tang, J., Yang, Z., Cao, Y., Ren, J., Benabu, Y., et al. (2016). Adaptive methylation regulation of p53 pathway in sympatric speciation of blind mole rats, Spalax. Proc. Natl. Acad. Sci. U.S.A. 113, 2146–2151. doi: 10.1073/pnas.1522658112
Keywords: sympatric speciation, environment heterogeneity, Ricotia lunaria, genetic divergence, reproductive isolation
Citation: Qian C, Yan X, Yin H, Fan X, Yin X, Sun P, Li Z, Nevo E and Ma X-F (2018) Transcriptomes Divergence of Ricotia lunaria Between the Two Micro-Climatic Divergent Slopes at “Evolution Canyon” I, Israel. Front. Genet. 9:506. doi: 10.3389/fgene.2018.00506
Received: 27 June 2018; Accepted: 08 October 2018;
Published: 14 November 2018.
Edited by:
Zhonghu Li, Northwest University, ChinaReviewed by:
Jing Wang, Norwegian University of Life Sciences, NorwayDongshi Wan, Lanzhou University, China
Copyright © 2018 Qian, Yan, Yin, Fan, Yin, Sun, Li, Nevo and Ma. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Eviatar Nevo, nevo@research.haifa.ac.il Xiao-Fei Ma, maxiaofei@lzb.ac.cn