Dynamic regulation of genome-wide pre-mRNA splicing and stress tolerance by the Sm-like protein LSm5 in Arabidopsis
© Cui et al.; licensee BioMed Central Ltd. 2014
Received: 29 June 2013
Accepted: 7 January 2014
Published: 7 January 2014
Sm-like proteins are highly conserved proteins that form the core of the U6 ribonucleoprotein and function in several mRNA metabolism processes, including pre-mRNA splicing. Despite their wide occurrence in all eukaryotes, little is known about the roles of Sm-like proteins in the regulation of splicing.
Here, through comprehensive transcriptome analyses, we demonstrate that depletion of the Arabidopsis supersensitive to abscisic acid and drought 1 gene (SAD1), which encodes Sm-like protein 5 (LSm5), promotes an inaccurate selection of splice sites that leads to a genome-wide increase in alternative splicing. In contrast, overexpression of SAD1 strengthens the precision of splice-site recognition and globally inhibits alternative splicing. Further, SAD1 modulates the splicing of stress-responsive genes, particularly under salt-stress conditions. Finally, we find that overexpression of SAD1 in Arabidopsis improves salt tolerance in transgenic plants, which correlates with an increase in splicing accuracy and efficiency for stress-responsive genes.
We conclude that SAD1 dynamically controls splicing efficiency and splice-site recognition in Arabidopsis, and propose that this may contribute to SAD1-mediated stress tolerance through the metabolism of transcripts expressed from stress-responsive genes. Our study not only provides novel insights into the function of Sm-like proteins in splicing, but also uncovers new means to improve splicing efficiency and to enhance stress tolerance in a higher eukaryote.
Immediately following transcription, many eukaryotic precursor messenger RNAs (pre-mRNA) are subjected to a series of modifications that are essential for the maturation, nuclear export and subsequent translation of these transcripts. One such modification, the removal (splicing) of non-protein-coding sequences from pre-mRNA, is an important step in gene regulation that also contributes to increased protein diversity from a limited number of genes. The precision and efficiency of splicing are critical for gene function . A non-precision splicing process would generate aberrant or non-functional mRNAs that are not only wasteful but can also lead to the production of unwanted or harmful proteins that may perturb normal cellular processes. Moreover, incorrectly spliced transcripts might also have a profound impact on other processes, including mRNA transcription, turnover, transport and translation. Accumulating evidence indicates that poor efficiency or defects in splicing can lead to diseases in humans [2, 3] and increase sensitivity to abiotic or biotic stresses in plants [4–6]. Although many molecular processes related to splicing have been well characterized, we still face a major challenge in understanding how precision and efficiency in splicing are regulated and how we could harness these regulations to enhance cellular functions.
Sm-like proteins (LSms) are a highly conserved family of proteins in eukaryotes both in terms of sequence and functions. LSms typically exist as heptameric complexes and play roles in multiple aspects of RNA metabolism [7–9]. The heptameric LSm1-7 cytoplasmic complex is located in discrete cytoplasmic structures called P-bodies, which are conserved in all eukaryotes and are thought to be involved in decapping and 5′ to 3′ RNA degradation [10, 11]. The LSm2-8 heptameric complex is located in the nucleus. This complex directly binds and stabilizes the 3′-terminal poly(U) tract of U6 small nuclear RNA, forms the core of the U6 small nuclear ribonucleoproteins (RNPs) and functions in pre-mRNA splicing [12, 13]. The Arabidopsis supersensitive to abscisic acid (ABA) and drought 1 (SAD1) gene locus encodes the LSm5 protein and was identified in a genetic screen for components that regulate the expression of stress-responsive genes in our previous work . SAD1 directly interacts with two other subunits, LSm6 and LSm7, and is a component of the LSm2-8 nuclear complex . Dysfunction of SAD1 increases the plant’s sensitivity to salt stress and to the stress hormone ABA in seed germination and root growth; moreover, sad1 mutants are defective in the positive feedback regulation of ABA biosynthesis genes by ABA and are impaired in drought stress induction of ABA biosynthesis, although the detailed molecular bases for these defects have not been identified. Recent studies suggested that the depletion of SAD1 and the other LSm protein (LSm8) reduced the stability of U6 RNPs and resulted in defects in pre-mRNA splicing that lead to intron retention in Arabidopsis [15, 16]. However, it is still unclear if the depletion of SAD1 or other LSm proteins has any effect on the selection of splice sites and alternative splicing (AS), and whether overexpression of these LSm proteins could affect splicing efficiency or accuracy.
To investigate possible regulatory roles of SAD1 protein in pre-mRNA splicing, we performed RNA sequencing (RNA-seq) of the wild-type Arabidopsis (C24 ecotype), the sad1 mutant and the SAD1-overexpressing plants (SAD1-OE). We found that SAD1 could dynamically control splicing efficiency and splice-site recognition and selection in Arabidopsis. Additionally, we discovered that SAD1 is required for regulation of splicing efficiency of many stress-responsive genes under stress conditions. Whereas there are increased splicing defects in sad1 mutants under salt-stress conditions, overexpression of SAD1 increases the splicing efficiency of stress-related genes. SAD1-OE plants are also more salt-tolerant than wild-type plants. Our work not only provides novel insights into the regulatory role of SAD1 and LSm proteins in splicing, but also suggests a new way to improve splicing efficiency and to optimize cellular functions and generate stress-resistant plants.
RNA sequencing of wild-type, sad1 mutant and SAD1-overexpressing plants
The Arabidopsis sad1 mutant was isolated in our previous genetic screen for components that regulate the expression of stress-responsive genes . The sad1 mutant was also more sensitive to stress and ABA inhibition of seed germination and seedling growth . Since loss-of-function mutations in any single-copied core LSm genes are expected to be lethal, the recovery of this point-mutation sad1 mutant provided an invaluable opportunity to study the functions of this important group of proteins. To explore the role of SAD1 in gene expression and stress tolerance, we generated transgenic Arabidopsis plants over-expressing the wild-type SAD1 cDNA (under control of the cauliflower mosaic virus 35S promoter) both in the wild type (ecotype C24) and in the sad1 mutant background (SAD1-OE, see Methods). Although the transgenic plants in both backgrounds had similar physiological and molecular phenotypes, here we mainly focus on SAD1-OE in the sad1 mutant background (referred to as SAD1-OE hereafter).
We next performed RNA-seq using the Illumina HiSeq platform (Illumina Inc., San Diego, CA, USA) on two-week-old seedlings of C24 (wild type), sad1 and SAD1-OE. These seedlings were subjected to two treatments: control (H2O) and salt stress (300 mM NaCl, 3 h). The salt-stress treatment was based on our previous observations that stress-responsive genes were most obviously activated and that the mutant sad1 showed strong molecular phenotypes under these conditions [14, 17]. Based on six cDNA libraries (C24-control, sad1-control, SAD1-OE-control, C24-NaCl, sad1-NaCl and SAD1-OE-NaCl), we generated a total of 164 million reads (101 bp in length, except for the SAD1-OE control, whose reads were 85 bp in length), about 90% of which could be uniquely aligned to the TAIR10 reference genome sequence (version TAIR10; ) (Additional file 1). Comparison of the mapped reads to the gene model (version TAIR10) revealed that approximately 95% of the reads mapped to the exonic regions, whereas only about 3% mapped to intergenic regions (Additional file 2), which were consistent with the Arabidopsis gene annotation. Plotting the coverage of reads along each transcript exhibited a uniform distribution with no obvious 3′/5′ bias, which reflects the high quality of the cDNA libraries (Additional file 3). Furthermore, assessing the sequencing saturation demonstrated that as more reads were obtained, the number of new discovered genes plateaued (Additional file 4). This suggests that extensive coverage was achieved, which can also be seen when the read coverage was plotted by chromosome, demonstrating extensive transcriptional activity in the genome (Additional file 5).
We previously identified the sad1 mutation as a G-to-A change 34 bp from the putative translation start site and predicted that the mutation would change a glutamic acid (E) residue to lysine (K). In the RNA-seq data, the mutation of sad1 at the genomic position 19,813,556 of chromosome 5 was confirmed. However, it turned out that the mutation occurred at the 3′ splicing acceptor recognition site of the first intron, changing the invariant AG dinucleotide to AA. Consequently, all of sad1 mRNAs were aberrantly spliced in the mutants, as visualized with the Integrative Genomics Viewer (IGV) browser [19, 20] (Figure 1C). We identified three main mutant transcripts in sad1: two with obvious aberrant 3′ splice sites (3′SSs) that respectively occurred 2 bp and 20 bp downstream of the mutated splice site; and one with the retention of the first intron (Figure 1C). All of these transcripts were validated by RT-PCR using primers that span the alternative 3′SSs, in which the corresponding events were detected in the sad1 mutant, but not in C24 (Figure 1D). Sequence analysis suggested that the transcript with aberrant 3′SSs that occurred 20 bp downstream of the mutated splice site did not alter the coding frame. It was predicted to produce one novel protein with the deletion of seven amino acids compared to the normal SAD1 protein. It seems that this mutant protein could provide some of the wild-type’s functions such that the sad1 mutation was not lethal. By contrast, the other alternative 3′SS and the intron retention led to a coding-frame shift that would generate a premature stop codon and thus would lead to truncated proteins. In the SAD1-OE plant, all these aberrantly spliced forms could be found, albeit at much lower levels than in sad1. However, normal SAD1 mRNA was overexpressed, with the transcript level more than 10-times higher than in C24, which was validated by quantitative RT-PCR (Figure 1E).
Identification of alternative splicing events in C24, sad1and SAD1-OE plants
To determine if there were any changes in pre-mRNA splicing upon the depletion or overexpression of SAD1, we first developed a pipeline to identify all AS events in C24, sad1 and SAD1-OE. The pipeline involved three steps: prediction of splice junctions, filtering of the false positive junctions and annotation of AS events. We randomly sampled 20 million uniquely mapped reads (estimated average approximately 57-times coverage on all the expressed transcripts) from each RNA-seq library for the identification or comparison of AS, respectively. This method ensured that the comparison of AS events would be performed at the same level.
To predict splice junctions, we mapped the RNA-seq reads onto the Arabidopsis genome using the software TopHat, which was designed to identify exon-exon splice junctions . After the alignment, we identified 732,808 junctions from the six RNA-seq libraries. Comparison of these junctions to the gene annotation (TAIR10) revealed that about 83% of total junctions had previously been annotated, and the remaining 17% were assigned as novel junctions (Additional file 6A). However, when trying to characterize these novel and annotated junctions, we found that there was a large number of novel junctions that had short overhangs (that is, fewer than 20 bp) with the corresponding exons, while most of the annotated junctions had large overhangs, with the enrichment at approximately 90 bp (Additional file 6B). Moreover, the novel junctions had relatively low coverage compared with the annotated junctions (Additional file 6C). In general, the junctions with short overhangs and lower coverage were considered as false positives, which are often caused by non-specific or error alignment. Therefore, to distinguish between true splice junctions and false positives, we assessed the criteria based on simulated data of a set of randomly constituted junctions. To do this, we first generated a set of 80,000 splice junctions in which annotated exons from different chromosomes were randomly selected and spliced together in silico. We also constructed 119,618 annotated junctions from the gene annotation. Since the length of our sequencing reads was 101/85 bp, the splice junction sequences were determined to be 180/148 bp long (90/74 nucleotides on either side of the splice junction) to ensure an 11 bp overhang of the read mapping from one side of the junction onto the other. Alignments to the random splice junctions were considered to be false positives, because such junctions are thought to rarely exist when compared to annotated junctions. The alignment of the raw RNA-seq reads to the random junctions revealed that 99.90% of false positive junctions had an overhang size of fewer than 20 bp (Additional file 7A). In sharp contrast, the alignment to the annotated junctions indicated that most (98.60%) annotated junctions had larger overhang sizes. In addition, we estimated that 56.90% of false positive junctions had only one read spanning the junction, whereas the annotated junctions had higher read coverage (Additional file 7B). To minimize the false positive rate, we required that the overhang size must be more than 20 bp and that there be at least two reads spanning the junctions. Using these criteria, we filtered out almost all the false positive junctions (Additional file 7C). Finally, we obtained a junction data set of 52,599 confident novel junctions from the six RNA-seq libraries. Based on these junctions, we identified all the AS events including cassette exons, alternative 5′SSs, alternative 3′SSs, mutually exclusive exons, coordinate cassette exons, alternative first exons and alternative last exons (Additional file 8).
Depletion of SAD1 activates alternative splicing
Twenty-two selected events were further validated by RT-PCR using the splicing-site-flanking primers, in which the corresponding AS events were detected in sad1 mutants, but were weakly or not presented in C24 (Figure 2C and Additional file 16). Figure 2C highlights three representative examples visualized by the IGV junction browser and validated by RT-PCR. The SBI1 (AT1G02100) gene had alternative 5′SSs in the 10th intron in sad1, but not in C24, an observation validated by RT-PCR using the forward primer that covered the splice junction and the reverse one that was located at the 11th exon. One can see that the corresponding isoform was detected in the sad1 mutant, but was not present in C24 (Figure 2C). The HINT3 (AT5G48545) gene had alternative 3′SSs in the fifth exon in the mutant sad1, which was validated by RT-PCR using a forward primer in the first exon and a reverse primer that covered the splice junction (Figure 2C). The gene PAC (AT2G48120) exhibited exon-skipping between the third and fifth exons, which was validated by RT-PCR using primers at the third and sixth exon, which meant that two different products were amplified, representing exon inclusion and skipping isoforms, respectively (Figure 2C).
Sequence analysis of these over-represented alternative 5′SSs and alternative 3′SSs (in the NaCl-treated sad1 mutant) revealed that these activated splice sites were still associated with GU and AG dinucleotides (Figure 2D; Additional file 17A), suggesting that the depletion of SAD1 did not change the accuracy of the sequence recognition of the splicing sites. When investigating the distribution of these activated splice sites, we found that alternative 5′SSs and 3′SSs were enriched in the downstream or upstream approximately 10 bp region of the dominant 5′SSs and 3′SSs, respectively (Figure 2E; Additional file 17B). This indicates that the depletion of SAD1 leads to the activation of the 5′SSs and 3′SSs proximal to the respective dominant ones. These results suggest that SAD1, as a component of U6 RNPs, may play a regulatory role in the selection of splice sites.
Interestingly, exon-skipping events also increased in sad1 mutants. When correlating each exon-skipping event with alternative 5′SSs and 3′SSs, we found that about 20% of the skipped exons simultaneously had alternative 5′SSs or 3′SSs in the mutants. This chance of occurrence was significantly higher than that expected for random sampling of all annotated exons (the probability of random occurrence was 0.02%, Fisher’s exact test, P <0.001). This result suggests a coordinated occurrence of exon-skipping and alternative splice site selection. Therefore, we considered that SAD1-depletion could simultaneously activate multiple alternative 5′SSs or 3′SSs that include not only the proximal ones, but also the distal ones, including those located at the next exons, albeit to a lesser extent. Nonetheless, the possibility that SAD1, probably as a splicing factor, may directly regulate exon-skipping in vivo could not be ruled out.
SAD1 depletion results in widespread intron retention
We next investigated if there is any influence of the splicing defects on the expression of the affected genes. Sequence analysis suggested that all of these intron retention events would generate premature stop codons in the intron-retained transcripts and, if translated, would produce truncated proteins. Although it is possible that some individual truncated proteins might still be functional, for our sequence analyses, we assumed that these intron-retained transcripts do not generate functional proteins. Through calculating the proportions of the intron-retained transcripts to the total transcripts for each gene with intron-retention in the mutant, we estimated that on average around 15% of total transcripts were with intron retention (Additional file 22). Moreover, when plotting the expression levels of the total and the functional transcripts (without intron) for each intron-retained gene between the wild type C24 and sad1 mutants (Additional files 23 and 24), we found that the expression levels of the total transcripts did not obviously change between C24 and sad1, but the functional transcripts tended to be down-regulated in the mutant. These results indicate that the splicing defects are associated with a global reduction of functional mRNAs, which might negatively affect the functions of these affected genes.
Genes with aberrant splicing in sad1are closely related to stress response and are activated by stress
Further analysis using Mapman  suggested that genes with aberrant splicing in sad1 mutants are involved in various stress response pathways, including hormone-signaling pathways, MAPK-signaling pathways and transcription regulation (Figure 4C; Additional file 29). Notably, some important genes (such as SnRK2.1 and 2.2, SOS2, DREB2A, NHX1, WRKY33, WRKY25, STT3A, CAX1 and RCI2A) involved in stress responses were identified to have splicing defects in the sad1 mutant. Among these genes, SnRK2.1 and 2.2 encode members of SNF1-related protein kinases activated by ionic (salt) and non-ionic (mannitol) osmotic stress that are required for osmotic stress tolerance ; SOS2 encodes a protein kinase essential for salt tolerance ; DREB2A encodes a transcription factor that activates drought and salt stress-responsive genes ; NHX1 encodes a vacuolar sodium/proton antiporter whose overexpression increases salt tolerance in several plant species including Arabidopsis; WRKY33 and WRKY25 encode plant WRKY transcription factors involved in response to salt and other stresses [30, 31]; STT3A encodes an oligosaccharyl transferase whose knockout mutants are hypersensitive to high salt conditions ; CAX1 encodes a high affinity vacuolar calcium antiporter and can be activated by SOS2 to integrate Ca2+ transport and salt tolerance ; and RCI2A (Rare-cold inducible 2A), whose product plays a role in preventing over-accumulation of excess Na+ and contributes to salt tolerance . These genes showed increased intron retention in the mutants, which were also validated by RT-PCR using intron-flanking primers where the corresponding intron-retained transcripts were more obviously identified in sad1, consistent with the RNA-seq data (Figure 4D). Above all, these results suggest that genes with aberrant splicing in sad1 are closely related to stress response, which could directly or indirectly contribute to the stress sensitivity of the sad1 mutant.
Overexpression of SAD1 rescues the splicing defects in the sad1mutant and strengthens splicing accuracy under salt stress
We further compared the expression levels of introns in SAD1-OE with those in C24 and sad1. We found that the expression of most introns in SAD1-OE was restored to normal levels (Figure 5F; Additional file 33), demonstrating that the intron retention indeed resulted from the sad1 mutation. Furthermore, using Fisher’s exact test we identified 76 introns from 75 genes that were significantly absent in NaCl-treated SAD1-OE, but were over-represented in NaCl-treated C24 (Additional file 34). This result shows that SAD1-overexpression can increase splicing efficiency.
Overexpression of SAD1improves plant salt tolerance
Further analysis suggested that these genes are involved in various stress response pathways (Additional file 37). Some of the stress-responsive genes that were more effectively spliced in SAD1-OE included ABF3/ABF2, encoding ABRE binding factors that mediate ABA-dependent stress responses [35, 36]; CIPK3, encoding CBL-interacting serine/threonine-protein kinase 3 that is involved in the resistance to abiotic stresses (for example, high salt, hyperosmotic stress) by regulating the expression of several stress-inducible genes ; and DREB2A, that encodes a transcriptional factor mediating high salinity- and dehydration-inducible transcription . These genes have been reported to be key regulators of ABA or salt-stress responses.
With the increased splicing efficiency in these key regulators of ABA or salt-stress responses, we were curious to know whether the SAD1-OE plants would have improved tolerance to salt stress. To test this, one-week-old seedlings of C24, sad1 and SAD1-OE grown on the regular Murashige and Skoog (MS) agar medium were transferred to MS agar plates supplemented with 0 (control), 50, 100 or 200 mM NaCl. We found that SAD1-OE seedlings showed enhanced tolerance to 100 mM NaCl on vertically placed plates (Figure 6C). At 200 mM NaCl, however, root elongation of all genotypes was inhibited and seedlings were not able to survive an extended period of the stress treatment (data not shown). Measuring the root growth of the seedlings showed that the roots of SAD1-OE were longer than those of C24 and sad1 at 100 mM NaCl (Additional file 38). We also tested salt tolerance of seedlings on horizontally placed agar medium plates. Two-week-old seedlings from ½ MS media were transferred onto 200 mM NaCl media and incubated for five days. The percentage of green leaf number over total leaf number was calculated for each seedling. The data indicated that SAD1-OE seedlings had a higher percentage of green leaves, suggesting that they were significantly less damaged by the salt stress than were the wild-type or sad1 seedlings (Figure 6D). To test further whether SAD1-OE plants were tolerant to salt stress at the adult stage and in soil, we grew these seedlings in soil and irrigated with either 50, 100, 150, 200 or 400 mM NaCl solutions at intervals of four days (see Methods). After two weeks of treatment, we found that sad1 plants were very sensitive to salt stress at concentrations above 150 mM and wild-type seedlings also exhibited signs of damages at higher salt concentrations as indicated by wilty inflorescence and damaged leaves, whereas the SAD1-OE plants were not obviously affected by the stress treatment and were also taller than the wild-type plants (Figure 6E; Additional file 39). These results indicate that SAD1-overexpression improves salt tolerance, which correlates with increased splicing accuracy and efficiency of stress-responsive genes.
Although studies in other eukaryotes, and more recently in plants, have demonstrated that LSm proteins 2-8, as the core of U6 RNPs, function in pre-mRNA splicing, whether or not these proteins have any roles in regulation of splicing efficiency and selection of splice sites has not yet been determined. In this study, through comprehensive transcriptome analysis of mutant and transgenic plants overexpressing the LSm5 gene SAD1, we demonstrated that SAD1 could dynamically regulate splicing efficiency and selection of splice sites in Arabidopsis. We also revealed that SAD1 modulates splicing of stress-responsive genes under salt-stress conditions. Finally, we showed that overexpression of SAD1 significantly improved splicing efficiency of the salt-responsive genes and resulted in enhanced salt tolerance in transgenic plants.
We found that SAD1 depletion activated the alternative 5′SSs and 3′SSs proximal to the dominant ones, suggesting that the wild-type SAD1 protein is necessary for precise splice-site recognition. To our surprise, relative to the wild-type plants, overexpression of SAD1 can strengthen the recognition accuracy and globally inhibit AS under salt-stress conditions. Therefore, we conclude that SAD1 can control selection of splice sites and splicing efficiency in a manner depending on SAD1’s abundance. This kind of splicing regulation, which could be referred to as a dynamic model, differs from but complements the kinetic model of splicing regulation [1, 38, 39]. In the kinetic model, the elongation rate of RNA polymerase II (Pol II) affects splicing efficiency such that a slower Pol II would allow more time for the recognition and processing of weak splicing sites so that splicing efficiency is enhanced. In the dynamic model, we reasoned that the spliceosome or other complexes involved in splicing are under thermodynamic equilibrium between association (complex formation) and disassociation (complex breakdown) at any given condition. A higher dosage of certain key small nuclear RNPs or splicing factors may drive the reaction toward the formation of the complex to enhance splicing efficiency. This dosage-dependent control of splicing suggests an alternative splicing regulation and it may be particularly important for the splicing of particular group of genes such as stress-inducible genes as discussed below.
Whereas increased alternative 5′SS usage was seen in sad1 both under control and salt-stress conditions (Figures 2A,B; Additional file 9A,B), the increase of alternative 3′SSs caused by SAD1 depletion and the inhibition of AS caused by SAD1 overexpression were only observed under salt-stress conditions. These findings indicate that SAD1 depletion or overexpression appears to impact splicing under salt-stress conditions more than under normal conditions. We considered that this distinct impact of SAD1 on splicing under normal versus stress conditions might have to do with increased transcriptional activation of the stress-responsive genes. Under salt-stress or other abiotic-stress conditions, plants activate the expression of a large number of stress-responsive genes that are not expressed or are expressed at lower levels under normal non-stressful conditions [40, 41]. With the simultaneous production of a large amount of these stress-inducible pre-mRNAs, cells would need to immediately recruit a significant amount of splicing factors and other factors for their co-transcriptional or post-transcriptional processing. This imposes a huge burden on the splicing machinery and as a result a significant portion of these transcripts fail to be processed adequately when the splicing machinery is compromised. This may be the reason why most of the splicing defective genes in sad1 are stress-regulated (Figure 4). Conversely, a higher SAD1 dosage could play a dominant role in enhancing the splicing efficiency of these salt-responsive genes through the promotion of recruitment and assembly of the splicing machinery as discussed above. As a result, the change in the AS pattern in SAD1-OE plants was more obvious under salt-stress conditions than under control conditions. Thus, the expression of these (and other) highly inducible genes may be particularly subjected to the dynamic regulation by certain splicing factors, which to some extent is similar to the kinetic regulation of splicing, both reflecting the saturated capability of cellular machinery.
We thought that the decreased splicing efficiency of the stress-responsive genes might contribute to the stress-sensitivity of the sad1 mutant. The splicing defects in sad1 lead to widespread intron retention in many stress-responsive genes (317 genes, Additional file 40). These genes include those encoding known key determinants of salt tolerance such as SnRK2.1/2.2, SOS2, DREB2A, NHX1, WRKY33, WRKY25, STT3A, CAX1 and RCI2A. The expression level of the functional transcripts for many of these genes were also found to be down-regulated (Additional file 41), although the cause of this down-regulation is unclear. All of these intron-containing transcripts were predicted to generate premature stop codons and truncated proteins if translated. This large-scale ‘hidden’ change in pre-mRNA splicing efficiency or gene expression, although relatively small for some of the individual genes, may collectively undermine plant’s readiness for the stress. However, it should be pointed that a direct relationship between the splicing defects and stress sensitivity in the sad1 mutant could not be established at this point.
Interestingly, an increase of splicing efficiency and expression of stress-responsive genes correlated with improved stress tolerance of the plants. Indeed, transgenic plants overexpressing SAD1 exhibited obviously increased tolerance to salt stress (Figure 6E), although the magnitude of the increase was moderate. Nonetheless, this finding is very significant for two reasons. First, it indicates that splicing efficiency may play an important role in regulating plant stress resistance. This is consistent with findings in several other genetic studies, where certain RNA processing factors were also found to be required for plant stress resistance. These factors include, for example, ABH1, LOS4 and RCF1, although the mechanisms involved were unclear. Secondly, our finding provides a new approach to improving plant stress resistance, namely, by regulating the splicing efficiency. Current methods to increase plant salt tolerance mostly involve the overexpression of structure genes such as ion transporter genes [29, 45, 46]. Constitutively expressing these structure genes may cause unwanted side effects that would result in reduced yield under normal growth conditions. However, enhancing splicing efficiency does not affect gene expression under normal conditions, as demonstrated in this study. Our finding may also be applicable to enhancing stress tolerance or other traits in other eukaryotic systems.
We demonstrated that SAD1 dynamically regulates splicing efficiency and plays a regulatory role in the selection of splice sites. Furthermore, we found that SAD1 specifically modulates splicing of the stress-responsive genes under stress conditions. Finally, we showed that overexpression of SAD1 improves salt tolerance of transgenic plant, which correlates with the increased splicing efficiency of the salt-stress-responsive genes. Our study provided novel insights into the regulatory role of SAD1 or LSm proteins in splicing and also suggested new strategies to improving splicing efficiency and bioengineering stress-resistant plants.
Materials and methods
Plant materials and growth conditions
The Arabidopsis sad1 mutant in the C24 background was described previously . For overexpressing the SAD1 gene, the SAD1 cDNA, amplified from the wild-type plant, was cloned into pENTRY1A. The LR reaction was then performed between pGWB502 and the pENTRY1A-SAD1. The resulting plasmid (pGWB502-SAD1) was introduced into Agrobacterium tumefaciens GV3101 and transformed into sad1 mutant plants using the floral dipping method. The transformants were selected on a ½ MS medium supplemented with 25 μg/ml hygromycin. Positive transformants were further confirmed by genotyping using the primers 5-CACCGGATCCTGATGGCGAACAATCCTTCACAGC-3; 5-TAATGAATTCGATCATTCTCCATCTTCGGGAGACC-3 for SAD1 cDNA. The confirmed transgenic seedlings (referred to as SAD1-OE) were used for RNA sequencing and RT-PCR analyses.
Seeds of C24, sad1 and SAD1-OE plants were sterilized with 50% bleach and 0.01% Triton X-100. The sterilized seeds were sown on ½ MS plates supplemented with 3% sucrose. After four-day stratification at 4°C, the plates were placed under a 16 h-light and 8 h-dark cycle at 21°C for germination and seedling growth. Twelve days later, the seedlings were treated with H2O (control) or 300 mM NaCl for 3 h, and harvested for total RNA extraction.
RNA extraction, library construction and sequencing
Using the TRIzol Reagent (15596–026, Invitrogen Co., Carlsbad, CA, USA), total RNAs were extracted from 12-day-old seedlings of wild-type C24, sad1 and SAD1-OE. Polyadenylated RNAs were isolated using the Oligotex mRNA Midi Kit (70042, Qiagen Inc., Valencia, CA, USA). The RNA-seq libraries were constructed using Illumina Whole Transcriptome Analysis Kit following the standard protocol (Illumina, HiSeq system) and sequenced on the HiSeq platform to generate high-quality single-end reads of 101 nucleotides (some with 85 nucleotides due to machine failure) in length.
RNA-sequencing data analysis pipeline
To analyze RNA-seq data, a pipeline was developed, which involved five steps: read alignment and junction prediction, the filter of false positive junctions, annotation of AS events, global comparison of AS and the identification of differential AS events (for details, see Additional file 42).
Read alignment and junction prediction
TopHat  was used to align the reads against the Arabidopsis genome sequences and annotated gene models were downloaded from TAIR10  allowing two nucleotide mismatches. Meanwhile, TopHat was also used to predict the splice junctions that did not permit any mismatches in the anchor region of a spliced alignment. The splice junctions were classified into known and novel splice junctions using the Perl script, which takes as input genome coordinates of all annotated exons and all predicted splice junctions. In addition, the expression levels of transcripts were measured by reads per kilo base per million values using the Cufflinks software .
The filter of false positive junctions
To estimate thresholds for filtering false positive junctions, two datasets of random and annotated splice junctions were first created. The dataset of 80,000 random splice junctions was created by joining each annotated 5′ donor sequence (90/74 bp from 5′SSs) and the annotated 3′ donor sequence (90/74 bp from 3′SSs) located on a different chromosome (Additional file 42). The 119,618 annotated splice junctions were created by joining each annotated 5′ donor sequence (90/74 bp from 5′SSs) and the annotated 3′ donor sequence (90/74 bp from 3′SSs) in order based on the gene annotation (Additional file 42). All splice junctions contained 90/74 nucleotides of exon sequence on either side of the junction to force an alignment overhang of at least 11 nucleotides from one side of the splice junction to the other. Then, the mapping software BWA  was used to align all reads to the random and annotated splice junctions that did not permit any mismatches. The alignments to random junctions were considered to be false positives, because such junctions are thought to rarely exist when compared to annotated junctions. We further characterized the false positive junctions, which generally have an overhang size of fewer than 20 bp and lower read coverage (Additional file 7A-B). To minimize the false positive rate, the overhang size with more than 20 bp and at least two reads spanning the junctions were required as cutoff value to filter the false positive junctions.
Annotation of alternative splicing events
JuncBASE  was used for annotating all AS events, including cassette exons, alternative 5′SSs, alternative 3′SSs, mutually exclusive exons, coordinate cassette exons, alternative first exons, alternative last exons and intron retention, which takes as input genome coordinates of all annotated exons and all confidently identified splice junctions. Notably, for identifying the events of intron retention, we required that at least five reads covered at least 50% of the region of one intron.
Global comparison of alternative splicing
The global comparison of AS among WT, sad1 and SAD1-OE was started with equally and randomly re-sampling uniquely-mapped reads to make sure that the comparison was at the same level. The comparison refers to the two facets: the absolute amount of each type of AS event and the number of junction reads that were assigned to each type of AS event, because both of them can be used to measure the global changes of AS. Meanwhile, Fisher’s exact tests in R  were used to identify differential representation of each type of AS event, performed on the number of junction reads that were assigned to each type of AS event.
The identification of differential alternative splicing events
Fisher’s exact tests were also used to identify differential representation of each AS event. For alternative 5′SSs and 3′SSs and exon-skipping events, Fisher’s exact tests were performed on the comparison of the junction-read counts and the corresponding exon-read counts between C24 and sad1 or SAD1-OE. The events with P <0.01 were identified as significantly different events. In addition, for those AS events that were uniquely identified in C24, sad1 or SAD1-OE, we would consider them significant if there were at least five junction reads to support and the P value of these events was assigned to equal zero. Similarly, for intron retention, Fisher’s exact tests were performed on the intron-read counts and the corresponding exon-read counts between C24 and sad1 or SAD1-OE. The events with P <0.001 were identified as significantly differential events. In addition, for those intron retention uniquely identified in C24 or the mutant, we would consider them significant if there was at least five-time coverage to support and the P value of these events was assigned to equal to zero.
The selected AS and intron retention events were validated by RT-PCR using a set of primers (Additional file 43) that were designed based on each AS event. Total RNAs from the C24, sad1 and SAD1-OE plants were extracted using Trizol solution (Invitrogen; cat.10837-08), treated with DNAase I, and reverse-transcribed to cDNA (random priming) by using a standard protocol (SuperScript II reverse-transcriptase, Invitrogen).
For the RT reaction, we used 3 μg total RNAs from the control (H2O) and 300 mM NaCl-treated C24, sad1 and SAD1-OE seedlings. The RT reactions were done with the Invitrogen SuperScript® III First-Strand Synthesis SuperMix in a 20 μl reaction system; the random Hexamer was used for first strand synthesis. The RT-solution was diluted 10 times, and 1 μl of the solution was used as template in 10 μl reaction system with 2 × SYBR Green (Invitrogen) Supermix (ROX). The quantitative RT-PCRs were performed in triplicate using the ABI 7900HT Fast Real-Time PCR System (Applied Biosystems Inc., Foster City, CA, USA). The primers 5-AAGGAGATAAG-GAGCTCGTTGG-3 and 5-ATCTGATCAAGCTTTGTGACC-3 were used for detecting expression levels of SAD1.
Salt-stress tolerance assays
Surface-sterilized seeds of C24, sad1 and SAD1-OE were sown onto agar plates with ½ MS and 1.2% agar. The plates were then kept at 4°C for four days before being incubated at 21°C for germination. Four days after germination, the seedlings were transferred to ½ MS agar plates supplemented with 0, 50, 100 or 200 mM NaCl, respectively. The seedlings were then allowed to grow for four days, and seedlings were photographed. Root length of these seedlings was measured by using Image J . For measuring leaf damage, two-week-old seedlings grown on ½ MS plates (0.6% agar) were transferred onto ½ MS medium plates supplemented with 200 mM NaCl and incubated for five days. The number of green leaves and yellowish or bleached leaves was counted for each seedling and percentage of green leaves among total leaves was calculated (leaves of the no-salt control treatment were all green and were not counted). To further test the tolerance to salt stress, seedlings grown on ½ MS agar plates were transferred to soil. One week after the transfer, the seedlings were irrigated with 50, 100 or 150 mM NaCl (in 1/8 MS salt), respectively . At 24 days after the transfer, the plants were further irrigated with 400 mM NaCl (100 ml) and pictures were taken four days later.
The RNA-seq data generated in this work has been submitted to the Sequence Read Archive database in NCBI. The accession number is SRP026082.
3′ splice sites
5′ splice sites
Murashige and Skoog
polymerase chain reaction
This work was supported by King Abdullah University of Science and Technology (KAUST), Kingdom of Saudi Arabia.
- Braunschweig U, Gueroussov S, Plocik AM, Graveley BR, Blencowe BJ: Dynamic integration of splicing within gene regulatory pathways. Cell. 2013, 152: 1252-1269. 10.1016/j.cell.2013.02.034.View ArticleGoogle Scholar
- Faustino NA, Cooper TA: Pre-mRNA splicing and human disease. Genes Dev. 2003, 17: 419-437. 10.1101/gad.1048803.View ArticleGoogle Scholar
- Singh RK, Cooper TA: Pre-mRNA splicing in disease and therapeutics. Trends Mol Med. 2012, 18: 472-482. 10.1016/j.molmed.2012.06.006.View ArticleGoogle Scholar
- Mazzucotelli A, Ribet C, Castan-Laurell I, Daviaud D, Guigne C, Langin D, Valet P: The transcriptional co-activator PGC-1alpha up regulates apelin in human and mouse adipocytes. Regul Pept. 2008, 150: 33-37. 10.1016/j.regpep.2008.04.003.View ArticleGoogle Scholar
- Palusa SG, Ali GS, Reddy AS: Alternative splicing of pre-mRNAs of Arabidopsis serine/arginine-rich proteins: regulation by hormones and stresses. Plant J. 2007, 49: 1091-1107. 10.1111/j.1365-313X.2006.03020.x.View ArticleGoogle Scholar
- Syed NH, Kalyna M, Marquez Y, Barta A, Brown JW: Alternative splicing in plants–coming of age. Trends Plant Sci. 2012, 17: 616-623. 10.1016/j.tplants.2012.06.001.View ArticleGoogle Scholar
- Beggs JD: LSm proteins and RNA processing. Biochem Soc Trans. 2005, 33: 433-438.View ArticleGoogle Scholar
- Collins BM, Harrop SJ, Kornfeld GD, Dawes IW, Curmi PM, Mabbutt BC: Crystal structure of a heptameric Sm-like protein complex from archaea: implications for the structure and evolution of snRNPs. J Mol Biol. 2001, 309: 915-923. 10.1006/jmbi.2001.4693.View ArticleGoogle Scholar
- Khusial P, Plaag R, Zieve GW: LSm proteins form heptameric rings that bind to RNA via repeating motifs. Trends Biochem Sci. 2005, 30: 522-528. 10.1016/j.tibs.2005.07.006.View ArticleGoogle Scholar
- Fenger-Gron M, Fillman C, Norrild B, Lykke-Andersen J: Multiple processing body factors and the ARE binding protein TTP activate mRNA decapping. Mol Cell. 2005, 20: 905-915. 10.1016/j.molcel.2005.10.031.View ArticleGoogle Scholar
- Ingelfinger D, Arndt-Jovin DJ, Luhrmann R, Achsel T: The human LSm1-7 proteins colocalize with the mRNA-degrading enzymes Dcp1/2 and Xrn1 in distinct cytoplasmic foci. RNA. 2002, 8: 1489-1501.Google Scholar
- Mayes AE, Verdone L, Legrain P, Beggs JD: Characterization of Sm-like proteins in yeast and their association with U6 snRNA. EMBO J. 1999, 18: 4321-4331. 10.1093/emboj/18.15.4321.View ArticleGoogle Scholar
- Verdone L, Galardi S, Page D, Beggs JD: LSm proteins promote regeneration of pre-mRNA splicing activity. Curr Biol. 2004, 14: 1487-1491. 10.1016/j.cub.2004.08.032.View ArticleGoogle Scholar
- Xiong L, Gong Z, Rock CD, Subramanian S, Guo Y, Xu W, Galbraith D, Zhu JK: Modulation of abscisic acid signal transduction and biosynthesis by an Sm-like protein in Arabidopsis. Dev Cell. 2001, 1: 771-781. 10.1016/S1534-5807(01)00087-9.View ArticleGoogle Scholar
- Perea-Resa C, Hernandez-Verdeja T, Lopez-Cobollo R, del Mar Castellano M, Salinas J: LSM proteins provide accurate splicing and decay of selected transcripts to ensure normal Arabidopsis development. Plant Cell. 2012, 24: 4930-4947. 10.1105/tpc.112.103697.View ArticleGoogle Scholar
- Golisz A, Sikorski PJ, Kruszka K, Kufel J: Arabidopsis thaliana LSM proteins function in mRNA splicing and degradation. Nucleic Acids Res. 2013, 41: 6232-6249. 10.1093/nar/gkt296.View ArticleGoogle Scholar
- Ishitani M, Xiong LM, Stevenson B, Zhu JK: Genetic analysis of osmotic and cold stress signal transduction in Arabidopsis: interactions and convergence of abscisic acid-dependent and abscisic acid-independent pathways. Plant Cell. 1997, 9: 1935-1949.View ArticleGoogle Scholar
- The Arabidopsis Information Resource. [http://www.arabidopsis.org]
- Thorvaldsdottir H, Robinson JT, Mesirov JP: Integrative Genomics Viewer (IGV): high-performance genomics data visualization and exploration. Brief Bioinform. 2012, 14: 178-192.View ArticleGoogle Scholar
- Robinson JT, Thorvaldsdottir H, Winckler W, Guttman M, Lander ES, Getz G, Mesirov JP: Integrative genomics viewer. Nat Biotechnol. 2011, 29: 24-26. 10.1038/nbt.1754.View ArticleGoogle Scholar
- Trapnell C, Pachter L, Salzberg SL: TopHat: discovering splice junctions with RNA-Seq. Bioinformatics. 2009, 25: 1105-1111. 10.1093/bioinformatics/btp120.View ArticleGoogle Scholar
- Huang DW, Sherman BT, Lempicki RA: Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat Protoc. 2009, 4: 44-57.View ArticleGoogle Scholar
- Huang DW, Sherman BT, Lempicki RA: Bioinformatics enrichment tools: paths toward the comprehensive functional analysis of large gene lists. Nucleic Acids Res. 2009, 37: 1-13. 10.1093/nar/gkn923.View ArticleGoogle Scholar
- Zimmermann P, Hirsch-Hoffmann M, Hennig L, Gruissem W: GENEVESTIGATOR. Arabidopsis microarray database and analysis toolbox. Plant Physiol. 2004, 136: 2621-2632. 10.1104/pp.104.046367.View ArticleGoogle Scholar
- Thimm O, Blasing O, Gibon Y, Nagel A, Meyer S, Kruger P, Selbig J, Muller LA, Rhee SY, Stitt M: MAPMAN: a user-driven tool to display genomics data sets onto diagrams of metabolic pathways and other biological processes. Plant J. 2004, 37: 914-939. 10.1111/j.1365-313X.2004.02016.x.View ArticleGoogle Scholar
- Fujii H, Verslues PE, Zhu JK: Arabidopsis decuple mutant reveals the importance of SnRK2 kinases in osmotic stress responses in vivo. Proc Natl Acad Sci USA. 2011, 108: 1717-1722. 10.1073/pnas.1018367108.View ArticleGoogle Scholar
- Liu JP, Ishitani M, Halfter U, Kim CS, Zhu JK: The Arabidopsis thaliana SOS2 gene encodes a protein kinase that is required for salt tolerance. Proc Natl Acad Sci USA. 2000, 97: 3730-3734. 10.1073/pnas.97.7.3730.View ArticleGoogle Scholar
- Sakuma Y, Maruyama K, Osakabe Y, Qin F, Seki M, Shinozaki K, Yamaguchi-Shinozaki K: Functional analysis of an Arabidopsis transcription factor, DREB2A, involved in drought-responsive gene expression. Plant Cell. 2006, 18: 1292-1309. 10.1105/tpc.105.035881.View ArticleGoogle Scholar
- Apse MP, Aharon GS, Snedden WA, Blumwald E: Salt tolerance conferred by overexpression of a vacuolar Na+/H+ antiport in Arabidopsis. Science. 1999, 285: 1256-1258. 10.1126/science.285.5431.1256.View ArticleGoogle Scholar
- Lai Z, Li Y, Wang F, Cheng Y, Fan B, Yu JQ, Chen Z: Arabidopsis sigma factor binding proteins are activators of the WRKY33 transcription factor in plant defense. Plant Cell. 2011, 23: 3824-3841. 10.1105/tpc.111.090571.View ArticleGoogle Scholar
- Jiang YQ, Deyholos MK: Functional characterization of Arabidopsis NaCl-inducible WRKY25 and WRKY33 transcription factors in abiotic stresses. Plant Mol Biol. 2009, 69: 91-105. 10.1007/s11103-008-9408-3.View ArticleGoogle Scholar
- Koiwa H, Li F, McCully MG, Mendoza I, Koizumi N, Manabe Y, Nakagawa Y, Zhu J, Rus A, Pardo JM, Bressan RA, Hasegawa PM: The STT3a subunit isoform of the Arabidopsis oligosaccharyltransferase controls adaptive responses to salt/osmotic stress. Plant Cell. 2003, 15: 2273-2284. 10.1105/tpc.013862.View ArticleGoogle Scholar
- Cheng NH, Pittman JK, Zhu JK, Hirschi KD: The protein kinase SOS2 activates the Arabidopsis H+/Ca2+ antiporter CAX1 to integrate calcium transport and salt tolerance. J Biol Chem. 2004, 279: 2922-2926.View ArticleGoogle Scholar
- Mitsuya S, Taniguchi M, Miyake H, Takabe T: Disruption of RCI2A leads to over-accumulation of Na+ and increased salt sensitivity in Arabidopsis thaliana plants. Planta. 2005, 222: 1001-1009. 10.1007/s00425-005-0043-9.View ArticleGoogle Scholar
- Kim S, Kang JY, Cho DI, Park JH, Kim SY: ABF2, an ABRE-binding bZIP factor, is an essential component of glucose signaling and its overexpression affects multiple stress tolerance. Plant J. 2004, 40: 75-87. 10.1111/j.1365-313X.2004.02192.x.View ArticleGoogle Scholar
- Yoshida T, Fujita Y, Sayama H, Kidokoro S, Maruyama K, Mizoi J, Shinozaki K, Yamaguchi-Shinozaki K: AREB1, AREB2, and ABF3 are master transcription factors that cooperatively regulate ABRE-dependent ABA signaling involved in drought stress tolerance and require ABA for full activation. Plant J. 2010, 61: 672-685. 10.1111/j.1365-313X.2009.04092.x.View ArticleGoogle Scholar
- Kim KN, Cheong YH, Grant JJ, Pandey GK, Luan S: CIPK3, a calcium sensor-associated protein kinase that regulates abscisic acid and cold signal transduction in Arabidopsis. Plant Cell. 2003, 15: 411-423. 10.1105/tpc.006858.View ArticleGoogle Scholar
- de la Mata M, Alonso CR, Kadener S, Fededa JP, Blaustein M, Pelisch F, Cramer P, Bentley D, Kornblihtt AR: A slow RNA polymerase II affects alternative splicing in vivo. Mol Cell. 2003, 12: 525-532. 10.1016/j.molcel.2003.08.001.View ArticleGoogle Scholar
- Kornblihtt AR: Coupling transcription and alternative splicing. Adv Exp Med Biol. 2007, 623: 175-189. 10.1007/978-0-387-77374-2_11.View ArticleGoogle Scholar
- Xiong LM, Schumaker KS, Zhu JK: Cell signaling during cold, drought, and salt stress. Plant Cell. 2002, 14: S165-S183. 10.1105/tpc.010278.View ArticleGoogle Scholar
- Yamaguchi-Shinozaki K, Shinozaki K: Transcriptional regulatory networks in cellular responses and tolerance to dehydration and cold stresses. Annu Rev Plant Physiol Plant Mol Biol. 2006, 57: 781-803.View ArticleGoogle Scholar
- Hugouvieux V, Kwak JM, Schroeder JI: An mRNA cap binding protein, ABH1, modulates early abscisic acid signal transduction in Arabidopsis. Cell. 2001, 106: 477-487. 10.1016/S0092-8674(01)00460-3.View ArticleGoogle Scholar
- Gong ZZ, Dong CH, Lee H, Zhu JH, Xiong LM, Gong DM, Stevenson B, Zhu JK: A DEAD box RNA helicase is essential for mRNA export and important for development and stress responses in Arabidopsis. Plant Cell. 2005, 17: 256-267. 10.1105/tpc.104.027557.View ArticleGoogle Scholar
- Guan QM, Wu JM, Zhang YY, Jiang CH, Liu RY, Chai CL, Zhu JH: A DEAD box RNA helicase is critical for pre-mRNA splicing, cold-responsive gene regulation, and cold tolerance in Arabidopsis. Plant Cell. 2013, 25: 342-356. 10.1105/tpc.112.108340.View ArticleGoogle Scholar
- Moller IS, Gilliham M, Jha D, Mayo GM, Roy SJ, Coates JC, Haseloff J, Tester M: Shoot Na+ exclusion and increased salinity tolerance engineered by cell type-specific alteration of Na+ transport in Arabidopsis. Plant Cell. 2009, 21: 2163-2178. 10.1105/tpc.108.064568.View ArticleGoogle Scholar
- Shi HZ, Lee BH, Wu SJ, Zhu JK: Overexpression of a plasma membrane Na+/H+ antiporter gene improves salt tolerance in Arabidopsis thaliana. Nature Biotechnol. 2003, 21: 81-85.View ArticleGoogle Scholar
- Trapnell C, Williams BA, Pertea G, Mortazavi A, Kwan G, van Baren MJ, Salzberg SL, Wold BJ, Pachter L: Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation. Nature Biotechnol. 2010, 28: 511-U174. 10.1038/nbt.1621.View ArticleGoogle Scholar
- Li H, Durbin R: Fast and accurate short read alignment with Burrows-Wheeler transform. Bioinformatics. 2009, 25: 1754-1760. 10.1093/bioinformatics/btp324.View ArticleGoogle Scholar
- Brooks AN, Yang L, Duff MO, Hansen KD, Park JW, Dudoit S, Brenner SE, Graveley BR: Conservation of an RNA regulatory map between Drosophila and mammals. Genome Res. 2011, 21: 193-202. 10.1101/gr.108662.110.View ArticleGoogle Scholar
- The R project for statistical computing. [http://www.r-project.org/]
- ImageJ. [http://rsbweb.nih.gov/ij/]
This article is published under license to BioMed Central Ltd. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.