1 Introduction

Photosynthesis is not only the most important chemical reaction to transform inorganic matter into organic matter through fixing solar energy on the earth but also the crucial source for plant carbon or biomass accumulation for growth and development. Its productivity is tightly linked to the delicate balance between carboxylation and oxygenation. As a bifunctional enzyme catalyzing both the carboxylation of D-ribulose-1, 5-bisphosphate (RuBP) that initiates photosynthetic CO₂ fixation and the oxygenation of RuBP that starts the nonessential photo-respiratory pathway (Nishimura et al., 2008), ribulose-1, 5-biphosphate carboxylase/oxygenase (Rubisco, EC4.1.1.39) plays a crucial role in the process of photosynthesis for most terrestrial plants. It consists of eight large subunit (LSU) and eight small subunits (SSU), compared to the little knowledge about the function of SSUs, LSUs have been proved to be the catalytic active site of the whole enzyme, which determines the catalytic efficiency; and the C-terminal of its amino acid was related to the fixation of CO₂.

Rubisco has been proved to evolve toward high-efficiency utilization of CO₂(Jordan and Ogren, 1981) for improving photosynthetic efficiency to adapt to environments (Galmes et al., 2005; Kubien et al., 2008). Previous studies showed that LSUs were encoded by the large single-copy region on the chloroplast genome, the rbcL gene (Strauss et al., 1988; Bausher et al., 2006). Although this gene was widely used as a molecular marker for the intergeneric or interfamilial level of phylogenetic analyses in angiosperms (Chase et al., 1993) with the characteristic of no evidence for adaptive evolution (Gould and Lewontin, 1979), as an important gene to encode the LSUs, rbcL was a likely target of natural selection to improve inefficient function of Rubisco, which brought the molecular evolution of the rbcL gene into question (Kapralov and Filatov, 2006). Recent researches showed that the evolution of the rbcL gene was indicated to be under positive selection in all principal lineages of land plants (Kapralov and Filatov, 2006). Furthermore, even under purifying selection in several species, such as Harveya purpurea, Striga gesnerioides, Orobanche fasciculata, and O. corymbosa (Wolfe and dePamphilis, 1997; Wolfe and dePamphilis, 1998; Leebens-Mack and DePamphilis, 2002). However, few studies focused on desert plants, which have faced extreme environment stresses such as the drought, UV, and high salt conditions for a long history.

Reaumuria soongarica (Tamaricaceae), a constructive xerophyte and halophyte shrub species, was widely distributed in all desert regions across the ACA (arid Central Asia) with an annual precipitation under 400 mm (Shi et al., 2013), including the Tengger Desert, Badain Jaran Desert, Gurbantunggut Desert, Gashun Gobi, Kumtag Desert, Qaidam Basin (on the northeastern Qinghai-Tibet Plateau, QTP), and Taklimakan Desert (which could be referred at www.eflora.org). Molecular phylogeographic study suggested that the evolutionary history of this Tertiary-relic species had been impacted by the uplift of the QTP and paleoclimate changes (Li et al., 2012; Yin et al., 2015). To be mentioned, as the typical C₃ plants (Ma et al., 2007), a recent trancriptome of R. soongoricapresented a complete set of C₄-function genes, which suggested that the regulation of Rubisco could be quite complicated for this plant to adapt to the extremely harsh environments. Furthermore, previous studies showed that the expression level and activase of the Rubisco were shifted in 142 land plants (Wang et al., 2011) under different environments, which also suggested a possible propensity to adapt to harsh habitats under different levels of salinity, water and temperature. Thus, it is reasonable to postulate that the rbcL subunit has undergone continual modification for better fitness in R. soongarica with shifts of paleogeography and paleoclimate.

In this study, to have a comprehensive acknowledge of whether the rbcL gene in R. soongarica experienced adaptive evolution, we first examined the sequence variations of the rbcL gene among 17 species of Tamaricaceae to seek the signal of adaptive evolution by applying phylogenetic analysis of maximum likelihood. Then, possible functional shifts were predicted by analysis of hydrophobicity and entropy on the locations under adaptive sites; the secondary and 3D-structure of Rubisco were also modelled to reveal higher levels of organization. Finally, to further clarify whether the expression level had shifted among different genotypes in R. soongorica in response to environment factors, the expression patterns of rbcL was analyzed under four treatments by real-time polymerase chain reaction (RT-PCR). Based on comprehensive analysis of the rbcL gene in R. soongorica, this study not only sheds light on the functional/structural features of Rubisco in R. soongorica but also provides useful information on directing genetic performance for enhancing the photosynthesis efficiency of desert plants to sustain fragile desert ecosystems, furthermore, to promote the ability to cope with desert aridification and global warming.

2 Material and methods 2.1 rbcL sequence acquirement and validation across biological database

The longest unigene with the annotation of ribulose-1, 5-biphosphate carboxylase/oxygenase was selected from the transcriptome of R. soongorica (Shi et al., 2013), and the peptide sequence was predicted on the GENSCAN Web Server at the Massachusetts Institute of Technology (MIT; http://genes.mit.edu/genscan.html). All the nucleotide sequences of rbcL from other species in Tamaricaceae were downloaded from the National Center for Biotechnology Information (NCBI; http://www.ncbi.nlm.nih.gov/). Finally, we obtained a total of 17 clean sequences (AY099899, AY099900, AY099902-AY099911, KC505169, KC505170, KC505172, Z 97650 and KM361003) in three genera, aligning with ClustalW implemented in Bioedit 7.2.2 (Hall, 1999). We also reconfirmed the ortholog gene of rbcL in R. soongorica by using the Basic Local Alignment Search Tool (BLAST) method on the Genebank database (http://blast.ncbi.nlm.nih.gov/Blast.cgi) and the Rubisco-large superfamily conserved domain by using Pfam (http://pfam.xfam.org/).

2.2 Phylogenetic analysis and adaptive selection detection

To confirm the phylogenetic relationships of rbcL in these 17 Tamaricaceae species, the maximum likelihood (ML) tree was reconstructed based on both the nucleotide and amino acid sequences analyses of rbcL using PhyML version 3.0 (Guindon and Gascuel, 2003) with 500 bootstraps under the best-fit model, respectively. The best-fit model for the nucleotide dataset was identified with the software of ModelGenerator, version 0.84 (Keane et al., 2004), based on Akaike information criterion (AIC) values (Akaike, 1974); and the TVM+I model and JTT+I model were finally chosen for nucleotide acid and amino acid sequences, respectively. Frankenia pulverulenta (HM850011) was used as the outgroup. Evolutionary changes and variations within and between species were judged with comparison of nonsynonymous-to-synonymous rates ratio (ω=dN/dS) by PAML package version 4.1 (Yang, 2007) with criteria as follows: when the ratio was equal (dN/dS=1), random drift of mutant alleles that are neutral at the molecular level was anticipated; when the ratio was lower than 1 (dN/dS < 1), purifying selection was expected; when the ratio was over 1 (dN/dS > 1), adaptive evolution or positive selection at the molecular level was indicated. The parameters of the M8 model using empirical Bayes approaches implemented in the CODEML program from PAML were used to calculate the posterior probabilities to identify amino acid sites potentially under positive selection.

2.3 Hydrophobicity and entropy analysis

To reveal the influence of mutation on biological diversity, and eventually contribute to the direction of biological diversity, the characteristics of the Shannon entropy of neighboring sites mutation were investigated with the program implemented in Bioedit 7.2.2 (Hall, 1999). Furthermore, to estimate the hydrophobicity force between polar AA side chains and the CO₂ penetration into Rubisco's active site, hydrophobicity was also analyzed with the program implemented in Bioedit 7.2.2 (Hall, 1999).

2.4 Structural analysis of Rubisco

As an alternative approach for functional annotation of novel protein sequences, Phyre2 algorithm (http://www.sbg.bio.ic.ac.uk/phyre2/) was used to analyze the structure of Rubisco. The confidence score of Phyre2 was established with 55% as the minimum cut-off value, and the proteins with confidence scores equal to or higher than this cut-off value are shown. The pdb file from phyre2 was visually displayed by PyMOL software (Schrodinger, 2010), and positive selection sites were also marked with the PyMOL software editor.

2.5 Total RNA Isolation and cDNA Synthesis and Quantitative Real-time PCR (RT-PCR)

To reveal the rbcL response to environmental factors, we germinated seeds of different geographic accession on Murashige-Skoog Medium (MS-medium) in the lab. After 14 days, the seedlings were exposed to drought, high-temperature, high-salt, and dark stresses. In detail, for drought stress, the seedlings were subjected to 15% PEG (polyethylene glycol) for 1 hour; for heat stress, the treatment was performed at 38 °C; for salt stress, plants were treated with 100 mmol/L NaCl for 3 hours; and for dark stress, plants were wrapped with foil paper for 2 hours. Total RNA of these materials was extracted using an E.Z.N.A. Plant RNA Kit (Omega Bio-tek, USA) with 2% PVBB added, and the remaining DNA was removed by RNase-free DNase (Omega Bio-tek, USA) according to the instruction manual. Total RNA concentration and purity was determined with ratio OD260/280 by NanoDrop1000 (Thermo Scientific), and all the samples passed quality control as the ratio of OD260/280 was between 1.9 to 2.2 and the ratio of OD260/230 < 2.0. RNA integrity was verified by 1.5% agarose gel electrophoresis with two clear bands of 28S/18S ribosomal RNA. For each sample, 1 μg of total RNA was reverse transcribed in a volume of 20 μL reaction with oligo dT primers by the RevertAidTM First Strand cDNA Synthesis Kit (Fermentas, Germany). The cDNAs were diluted by 1:10 with nuclease-free water prior to the qPCR analyses.

RT-PCR primers for each unigene sequences fromde novo transcriptome ofTamarix hispida were designed by OligoPerfect Designer (http://tools.invitrogen.com/content.cfm?pageid=9716&icid=fr-oligo-6), with the length approximately 20 bp, optimal melting temperatures 60 °C, and optimal GC content 50%. Finally, two pairs of primers (rbcL and H2A) from Yan et al. (2014) with a product size between 100 to 300 bp were used in this study. PCR efficiency for each sample was calculated by LinRegPCR, with the mean PCR efficiency per amplicon between 1.8 and 2.0. A total of 20 μL reaction-system volume was used for amplification, consisting of 10 μL DyNAmo Flash SYBR Green qPCR Kit Master Mix (Thermo Scientific), 0.5 μL of forward and reverse primer, respectively, 0.2 μL of F-402 buffer, 2 μL cDNA synthesized from total RNA, and 6.8 μL double-distilled water. The PCR program contained an initial denaturation step of 5 min at 95 °C, followed by denaturation for 15 s at 95 °C, annealing for 30 s at 60 °C, and extension for 30 s at 72 °C for 40 cycles. The real-time PCR thermal cycler qTOWER 2.0/2.2 (Analytik Jena, Germany) was used to obtain relative expression levels of each sample. The dissociation curve was obtained by heating the amplicon from 60 °C to 95 °C for each primer to get good specificity and quality, efficiency, and the dissociation curve. The Cq value per sample and the fluorescence threshold was set as determined by LinRegPCR, and relative quantification was determined using the Delta Delta Ct method. The relative transcripts were normalized with reference genes from Yan et al. (2014).

3 Results 3.1 Gene identification and sequence analysis

The sequence of Unigene19142_A from the transcriptome of R. soongorica(Shi et al., 2013) was annotated as ribulose-1, 5-bisphosphate carboxylase/oxygenase large subunit. The predicted peptide on GENSCAN was 475 amino acids with 6 bp polyA at the position -289 to -284 and 5'-UTR at the position -2153 to -2114, which was similar to the previously published sequence of R. soongorica (AAM26908) with 100% identity, except for the decrease of 25 amino acids at the N-terminal and 28 amino acids at the C-terminal. Pfam results proved that the Rubisco-large superfamily was conserved even in different families (Figure 1).

3.2 rbcL phylogenetic tree of three genera in Tamaricaceae

A total of 17 clean rbcL sequences with 1, 200 bp and their predicted protein sequences from Tamaricaceae, were used to reconstruct the nucleotide and protein phylogeny trees (Figure 2). As shown in Figure 2a, the tree is based on the variation of DNA sequences: the three genera, Reaumuria, Myricaria, and Tamarix, were independently separated from the main tree. In the positions of 336 and 595 of the DNA fragment, the mutations "T" and "G" were fixed in the Tamarix and Reaumuria genera, respectively, whereas, mutation "G" at position 660 was still in a steady-state fitness variation. However, the topology of the protein phylogeny tree was much different, as shown in Figure 2b; although the value of the bootstrap was low, Reaumuria was separated from the clade of Myricaria and Tamarix, and the species of these two genera were complicated to get a relatively clear lineage relationship. The discrepancy between the nucleotide and protein phylogeny trees suggested the level of differentiation in Tamaricaceae was higher in the transcript and translation process, which may contribute more to the divergence of Tamaricaceae species.

3.3 Characteristics of selective sites of rbcL in R. soongorica

Results of the nonsynonymous/synonymous rate ratio ω showed that, with posterior probability over 95% by empirical Bayes analysis, three mutations at positions 112, 199, and 220 on rbcL were under positive selection (marked with an asterisk in Table 1). At position 112, the amino acid alanine was substituted by serine (Ala¹¹²→Ser¹¹²), which was largely harbored in the genera Tamarix and Reaumuria; at position 199, the amino acid threonine was replaced by alanine or glycine (Thr¹⁹⁹→Ala¹⁹⁹ or Thr¹⁹⁹→ Gly¹⁹⁹), which was largely harbored in genera Reaumuria and Myricaria; whereas, at position 220, the isoleucine was replaced by methionine (Ile²²⁰→ Met²²⁰), which was harbored in all the three genera (Table 2). As shown in Figure 1, all these substitutions were found in the terminal groups rather than the base group. According to the nucleotide sequences, the fixed "T" at position 336 was related to the amino acid mutation of Ala¹¹²→Ser¹¹², the mutation "G" at 595 was related to the mutation of Thr¹⁹⁹→Ala¹⁹⁹ or Thr¹⁹⁹→Gly¹⁹⁹, while the complex mutation at 660 loci was clearly evolved from isoleucine to methionine under positive selection. Furthermore, the composite structure analysis showed that the sequence of R. soongorica harbored 11 substitutions, including eight non-selective positions and three selective positions (Figure 3a), and the 3D structure for Rubisco in R. soongorica showed that all of the three positive-selection positions (112, 199, and 220) were located in the alph-helix of the α/β barrel domain (red marker in Figure 3b).

3.4 Hydrophobicity and entropy analysis of selective site in R. soongorica

From the Eisenberg scale, mean hydrophobicity analysis by Bioedit 7.2.2, the hydrophobicity decreased and polarity enhanced at the three selective loci (Figure 4a), suggesting that the hydrophobicity force between the polar AA side chains was shifted; and the conformation of the protein tends to be changed. Furthermore, the entropy was also increased, ranging from 0.6 to 1 (Figure 4b), also suggesting that those sites were unstable. The decreased hydrophobicity and increased entropy presumably facilitate CO₂ penetration into Rubisco's active site, which suggests those selection loci were tightly impacted by a great change of CO₂ concentration with the paleoclimate shift, to promote the evolution at the expense of a considerable energy cost to adapt to the environment.

3.5 Response of rbcL to environmental factors

Results of the RT-PCR showed that the average expression level of rbcL in the HG accessions was slightly upregulated under various stresses, while in the XGG accessions they were somewhat downregulated (Figure 5). Under the PEG treatment, the expression levels of rbcL were slightly increased and were higher in XGG accessions than HG ones. Under the heat treatment, the expression level of rbcL was slightly depressed or showed no change in HG accessions, whereas it was slightly increased in the accession of XGG. Under the dark treatment, the transcript level of rbcL remained unchanged except XGG030. Carbon assimilation is one of the major sites sensitive to drought, heat stress, and light in photosynthesis. The slight regulation of rbcL expression under PEG, heat, and dark treatments revealed that the maintenance mechanism of CO₂-assimilation capacity was mainly protected by Rubisco activity or degradation of Rubisco. Different from the other three accessions, the HG020 elevated to a high expression level under the salt treatment. This finding could be attributed to both the genetic variation in the local environment and the regulation of rbcL expression.

4 Discussion and conclusion

Although, in the field of phylogenetic analyses, the rbcL gene has been considered to be conserved (Chase et al., 1993), issues of rbcL evolution under natural selection had been proposed recently within all principal lineages of land plants. Most of the Rubisco residues appeared to be under purifying selection (Wolfe and dePamphilis, 1997; Wolfe and dePamphilis, 1998; Leebens-Mack and DePamphilis, 2002), while only a small fraction have been suggested to be under positive selection in particular taxonomic groups (Kapralov and Filatov, 2006; Christin et al., 2008; Iida et al., 2009; Miwa et al., 2009; Sen et al., 2011; Wang et al., 2011; Young et al., 2012). These phenomena indicated that the evolution of rbcL, which was induced by the constant fine-tuning of Rubisco to adapt to different environments, was always occurring (Galmés et al., 2005).

In this study, three nonsynonymous mutations were suggested to be under positive selection in the family Tamaricaceae (Table 1). They were harbored in the terminal groups rather than the base group (Figure 1); and the fixed mutation of "T" at position 336 contributed to the Ala¹¹²→Ser¹¹²; the mutation "G" at position 595 contributed to Thr¹⁹⁹→Ala¹⁹⁹ or Thr¹⁹⁹→ Gly¹⁹⁹; and the flexible mutation at position 220 was evolved from isoleucine to methionine under positive selection. On the level of nucleotide acid, the mutations "T" and "G" were fixed in Tamarix and Reaumuria genera, respectively, whereas, mutation "G" at loci 660 was still in a steady-state fitness variation, which contributed to the independent clade for the three genera, Reaumuria, Myricaria, and Tamarix. On the level of amino acid, Ala¹¹²→Ser¹¹²was harbored majorly in genera Tamarix and Reaumuria; Thr¹⁹⁹→ Ala¹⁹⁹ or Thr¹⁹⁹→Gly¹⁹⁹ was harbored majorly in genera Reaumuria and Myricaria, whereas, Ile²²⁰→ Met²²⁰ was harbored in all three genera. In considering their habitat environments, we suggest that these replacements may account for the ecological differences; and Ala¹¹²→Ser¹¹²seemed to be more responsible for the drought-and salinity-selection pressure, as water is demanded for increasing survival and the salt-secreting ability was reduced from Reaumuria to Tamarix toMyricaria. The spatial analysis showed that all the three positive-selection loci (112, 199, and 220) were located in the alph-helix (red marker in Figure 3b), and previous studies showed that the positively selected sites preferably located in the helices (Wang et al., 2011). Furthermore, we also found that these substitutions were introduced to the decrease of hydrophobicities and the increase of the entropy (Figure 4), supporting the contention that these selective substitutions could induce the instability of the Rubisco protein structure and contribute to the diversity of species. According to Yin et al. (2015), the intrafamily divergence of Tamaraceae was dated to approximately 66.36 Ma, and the two genera Tamarix and Myricaria could have diverged at approximately 43.26 Ma, when a higher concentration of CO₂ was in the atmosphere. The high concentration of CO₂ induced the faster evolution of Rubisco to more efficient utilization of CO₂ (Jordan and Ogren, 1981; Savir et al., 2010). As an example, R. soongarica had been regarded as the Tertiary-relic species (Hou, 1987), which represented the background of a great shift in global CO₂ concentration (Petit et al., 1999; Zachos et al., 2001) around Rubisco. At that time, C₄ plants arose from C₃ plants via the gradual addition of constituents and their ability to thrive across a diversity of habitats (Christin and Osborne, 2013). The variation of δ¹³C value for R. soongorica ranged from -22.8‰ to -29.9‰, and the lowest δ¹³C appeared to favor salinity. What's more, the values increase when the salinity is lower or higher than the optimum level; and the variation of δ¹³C may be attributed to all of the genes encoding key enzymes in the C₄carbon-fixation pathway, detected in the transcriptomic dataset by Shi et al. (2013); thus, we would like to infer that there may exist the intermediate of C₃-C₄ photosynthesis in R. soongorica to adapt to different environments, as found in other species (Ueno, 1998; Brautigam et al., 2011; Williams et al., 2012).

Environmental factors such as salinity, temperature, and drought could induce the limitation of CO₂ fixation and further result in the decrease of photosynthetic capacity in plants (Barhoumi et al., 2007; Hozain et al., 2010; Carmo-Silva and Salvucci, 2012; Du et al., 2015; Guo et al., 2015), which were suggested to be a driving force for the evolution of Rubiscio (Jordan and Ogren, 1981; Delgado et al., 1995; Galmes et al., 2015). As found in other plants such as the soybean (Carmo-Silva et al., 2010) and tomato (Wang et al., 2015), the transcriptional levels of rbcL in R. soongorica, which could reflect the maintenance of CO₂-assimilation capacity for plants, was slightly regulated by dark, heat, and PEG treatments (Figure 5). However, the response to the salt treatment was complex, as shown in Figure 5, with the expression level of rbcL highly upregulated in HG020 accession, while no change or slight deregulation in the other three accessions. In fact, the expression level of rbcL is not only associated with the concentration of Rubisco (Randle and Wolfe, 2005; Mehrotra et al., 2011) but also with the activation state of the enzyme (Galmes et al., 2013) and the synthesis and degradation of Rubisco or the regulation mechanism (Mori et al., 2002; Bowman et al., 2013; Wang et al., 2015). Until now, the precise modulation of Rubisco activity under varying environmental conditions was unclear, such as different responses to drought stress were found in different species (Sharkey and Seemann, 1989; Parry et al., 2002; Galmes et al., 2013; Zhou et al., 2015).

Numerous phylogeographic studies have suggested that natural selection would drive adaptive diversification of genes (Liu and Zhu, 2010; Zhao et al., 2014; Huang et al., 2015; Sun et al., 2015) and even result in intraspecific variation (Fu et al., 2013; Johnson et al., 2015; Yin et al., 2015). A great genetic variability in salt tolerance among genotypes of species had been shown (Lowry et al., 2009; Bchini et al., 2010), and salinity tolerance was positively related to the genetic diversity of species (Zhang and Chen, 2010). In this study, we found that the expression level of the Rubisco and its activase were shifted in response to variation of CO₂ concentration in the history; furthermore, the heterogeneous desert environments may also contribute to the functional divergence of Rubisco. To have a deeper understanding of the adaptive evolution of rbcL, further experiments with wider sampling are still needed. At present, this study presented an example of the way that Tamaricaceae plants orchestrate the C₃ and C₄ pathways to adapt to environmental changes, and it further sheds light on the evolution scenarios from the C₃-like pathway to the C₄ pathway.

Acknowledgments:

This work was supported by the National Natural Science Foundation of China (NSFC, Grant Nos. 31370395 and 31500266) and the "One Hundred Talents" project of the Chinese Academy of Sciences (Grant No. 29Y127E71).