Alternative splicing produces different mature RNAs from a single gene. Its impact on increasing transcript diversity has continued to broaden our understanding of gene regulatory mechanisms since it was first observed in 1977 1,2 . The potential to create novel transcript diversity via alternative splicing is immense. The Drosophila DSCAM gene, which functions as an axon guidance receptor, is an extreme example of alternative splicing. It contains 115 exons and is estimated to give rise to more than 38,000 isoforms that are spatiotemporally regulated to achieve specific regulation in Drosophila neural development 3,4 . Alternative splicing events can be observed in more than 95% of human genes 5 . High throughput proteomics and ribosome bound mRNA sequencing (Ribo-Seq) studies show that a significant fraction of alternative splice variants are translated into protein isoforms 6,7 . Additionally, an ever-increasing array of transcriptome sequencing has revealed the existence of novel non-coding RNAs generated through alternative splicing suggesting their importance in gene regulatory circuits in all eukaryotic clades 8,9 . Differential splicing in closely related species have shown to reflect their divergent adaptive strategies not readily detectable at the primary gene expression level 10 . Tissue-and species-specific splicing is more divergent than promoter level divergence among closely related species facilitating independent evolutionary trajectories fitting to each species as highlighted by Calarco et al. 11 Therefore, genome wide discovery of new transcripts produced via alternative splicing becomes a critical initiative to understand the diversity of gene products and systematically assess their role in evolutionary innovations.
(which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for this preprint this version posted October 17, 2022. ; https://doi.org/10.1101/2022.10. 13.512046 doi: bioRxiv preprint 3
Alternative splicing increases proteome diversity as well as regulatory complexity in plants [12][13][14] . In the model plant Arabidopsis thaliana, majority of the genes (>60%) undergo alternative splicing. There are more than 70,000 non-redundant transcript isoforms reported for A. thaliana 15,16 . Similar reports on maize 17 , sorghum 18 , and cotton 19 demonstrate that alternative splicing is prevalent in plants. Differential splicing has also been targeted in crop breeding as shown with sunflowers 20 .
Large scale changes in alternative splicing have been reported to allow transcriptional adjustments in response to abiotic stresses including salt 21 , cold 22 , hypoxia 23 , and heat stress 24 . Targeted functional studies have also highlighted the significance of alternative splicing in responses to abiotic stresses. For example, the heat shock protein gene, hsf2 in A. thaliana produces an alternatively spliced transcript resulting in a truncated protein that in turns binds to the hsf2 promoter to enhance transcription of hsf2 during heat stress 25 . While the majority of published studies converge on alternative splicing being a key mechanism for environmental stress adaptation in plants, such studies are limited to abiotic stress sensitive crop models or to A. thaliana.
Compared to crop plants, extremophytes that are naturally found in extreme environments are equipped with evolutionary innovations that give them the ability to cope with multiple and extreme levels of environmental stresses 26 . Therefore, extremophytes could show how the expanded transcriptome diversity via alternative splicing may render additional paths for abiotic stress adaptations absent in stress sensitive models.
In this study, we have used the extremophyte model, Schrenkiella parvula (formerly Thellungiella parvula and Eutrema parvulum) 27,28 . to examine its transcriptome diversity augmented by alternative splice variants. S. parvula shares a highly co-linear genome with A. thaliana 29,30 . Yet, S. parvula is uniquely adapted to multiple abiotic stresses reflecting its natural habitats often associated with hypersaline lakes in the Irano-Turanian region 31,32 .
Previous studies have shown that the S. parvula genome is enriched with duplicated genes associated with abiotic stress responses and stress responsive genes show constitutive high expression as a stress preadaptation compared to A. thaliana 29,32,33 . Alternative splicing plays a complementary role to gene duplications and provides an additional path to increase transcript diversity 34 . Therefore, we aimed to test the overall hypothesis that alternative splicing leads to increased diversity of stress responsive transcripts in S.
In this study we investigated the complexity of the alternative splicing landscape in roots and shoots in response to salt stress and how alternative splicing may provide transcript diversity associated with adaptations to environmental stress in the model extremophyte S. parvula. We used PacBio Iso-Seq sequencing to identify and annotate alternative splice variants and Illumina short reads to quantify isoform abundance. We find that the S. parvula transcriptome is enriched in stress-associated isoforms. It shows specific isoform usage in a less disordered state compared to the stress-sensitive model A. thaliana in response to high salinity.
(which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for this preprint this version posted October 17, 2022. ; https://doi.org/10.1101/2022.10.13.512046 doi: bioRxiv preprint 4 Schrenkiella parvula (ecotype Lake Tuz, Turkey; Arabidopsis Biological Resource 575 Center/ABRC germplasm CS22663) seeds were grown hydroponically as previously described 33 . Briefly, plants were grown at a light/dark cycle of 12/12 hr, 100 -120 mM•m -2 s -1 photon intensity, 20-22 °C, and 1/5 th Hoagland's solution for four weeks. These were treated with a with a combination of 250 mM NaCl, 250 mM KCl, 30 mM LiCl, and 15 mM H 3 BO 3 for three days to generate tissue samples used to create a reference transcriptome with PacBio Iso-seq sequencing. Shoots and roots were harvested separately. RNA was extracted using QIAGEN RNeasy Plant Mini Kit (QIAGEN, Hilden, Germany) with column digestion to remove DNA contamination. About 4 µg of total RNA per tissue type at a quality of RNA integrity number ≥ 8 based on a Agilent 2100 Bioanalyzer (Agilent Technologies, CA, USA) were used to generate RNA-Seq libraries.
Shoot and root (1 μg) extracted as described above were used for cDNA synthesis using the SuperScript cDNA Synthesis Kit (Invitrogen, Massachusetts, USA) following manufacturer's instructions to test the presence of multiple isoforms independent from Iso-Seq for a randomly selected gene set expected to express multiple isoforms. Isoform specific PCR primers (Supplementary Table 1) that span the alternative splice sites were designed to use with an amplification protocol (initial denaturation at 95 °C for 3 min; 30 cycles of 95 °C for 30 s, 50-56 °C for 30 s, 72 °C for 2.30 to 3 min ; 72 °C for 10 min) run on a Bio-Rad T100 Thermal Cycler (Hercules, CA, USA) with a PCR Master mix Solution i-MAX II (iNtRON Biotechnology, S Korea). PCR products were separated on a 1% agarose gel.
To quantify isoform abundance in response to high salinity compared to control conditions, RNA was extracted from hydroponically grown S. parvula and A. thaliana (Col-0) as described in Tran et al. (2021).
These plants were treated with 150 mM NaCl for 24 hours and harvested together with samples hydroponically grown without added NaCl as a control condition. The hydroponic growth conditions except for the specific salt treatment was kept equivalent to growth conditions given to plants used for reference transcript generation with Iso-Seq. At least 5 plants were used per biological replicate and three biological replicates were used for each root and shoot sample for S. parvula and A. thaliana to yield a minimum of 1 µg of total RNA per sample used for standard RNA-seq library preparation.
For Iso-Seq based long read sequencing, cDNA synthesis, sequencing library preparation, and PacBio sequencing were conducted at the Arizona Genomics Institute, University of Arizona, USA. Two Iso-seq sequencing SMRT libraries were constructed following size selection from ≤ 4 kb and ≥ 3.5 kb per each tissue and ran on two Pacific Biosciences Sequel cells with v2.1 Chemistry. For RNA-seq based short read sequencing, mRNA enriched cDNA synthesis, library preparation, and sequencing were conducted at the Roy K. Carver Biotechnology Center, University of Illinois Urbana-Champaign, USA. Briefly, True-Seq strand specific libraries (Illumina, San Diego, CA, USA) were multiplexed and sequenced on an Illumina HiSeq4000 platform to generate >15 million 50-nucleotide single-end reads per sample.
(which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for this preprint this version posted October
Raw Sequel data were processed using isoseq_sa5.1 pipeline (https://github.com/PacificBiosciences/IsoSeq_SA3nUP. Circular consensus sequences (CCS) were generated from subread BAM files with following parameters: minLength = 50, -noPolish --minLength=50, --maxLength=15000, --minPasses=1, --minPredictedAccuracy=0.8, --minZScore=-999 --maxDropFraction=0.8. CCS reads were selected as full length reads if it contained the 5′ and 3′ primers and a poly(A) signal preceding the 3′ primer without additional copies of adapters. The full length consensus transcripts were further clustered using ICE (Iterative Clustering for Error Correction) to obtain high-quality isoforms with post-correction accuracy above 99% using Quiver. Error corrected full length reads were mapped to the Schrenkiella parvula reference genome v2.2 (Phytozome genome ID: 574) to annotate isoforms assigned to gene models and further select a set of high confidence transcript models. An isoform is annotated as a full length transcript mapped to a genomic locus that has a single gene model assigned. If more than one isoform is mapped to a gene model, the second and subsequent isoforms are considered products of alternative splicing.
First, TAPIS 17 was used to map isoforms and further error correct the isoforms. To map reads to the genome GMAP 36 was used with parameters, --no-chimeras, --cross-species --expand-offsets 1, -K 3000. Then, SQANTI 37 was used with default parameters to identify the isoform that matched the primary gene model in the genome and to assign additional isoforms that may be derived from that gene model as alternatively spliced isoforms if both types of full-length isoforms were present in our processed full length data. Canonical splice sites were defined as AG at the acceptor site and GT at the donor site. All the other splice sites were categorized as non-canonical splice sites. Custom python script was used to count canonical and non-canonical splice sites. Finally, we selected non-redundant structurally distinct isoform models that also contained a complete and uninterrupted open reading frame as a selected set of putative protein coding transcript models.
Only isoforms that are likely to code for proteins were used for downstream analyses in the current study due to the high uncertainty of functional significance and limited annotation resources available for newly identified non-coding isoforms.
Functional annotations were assigned using PANTHER 38 and A. thaliana Gene Ontology (GO) annotations (version release date 2020-07-16; DOI:10.5281/zenodo.3954044). Test for enriched functions were performed using BiNGO 39 . Further clustering of enriched functions were performed using GOMCL 40 (with parameters: -gosize 1500 -Ct 0.7 -I 1.5 -hm -nw -d -hg 0 -hgt -ssd) to get a non-redundant set of representative functional annotations at p-values ≤0.05 adjusted for false discovery rate.
Following quality checks using FastQC (http://www.bioinformatics.babraham.ac. uk/projects/fastqc/). RNA-seq reads were mapped to gene models for A. thaliana (TAIR10) or S. parvula (Reference v2.2) as well as transcript models obtained from AtRTD2 41 or S. parvula Iso-seq supplemented transcript models using pairs between S. parvula and A. thaliana were assigned based on Oh & Dassanayake 2019 43 . RNA-seq reads mapped to gene models were used to identify differently expressed genes. A custom python script was used to count uniquely mapped reads to each gene model. Differentially expressed genes between control and salt treatments within each species were identified using DESeq2 44 RNA-seq reads mapped to transcript models were used for generating expression values for isoforms as well as quantify alternative splicing event frequency. Expression counts for isoforms were converted to TPM (Transcript Per Million) and in comparisons where an isoform was counted as expressed had ≥ 0.5 TPM normalized expression per isoform independent from the expression quantified at the gene level. Isoform ratio per ortholog pairs was calculated based on the number of isoforms per S. parvula ortholog divided the number of isoforms detected in the A. thaliana ortholog.
Differential splicing was assessed using SUPPA2 45 . Briefly, alternative splice events were identified using generateEvents program and differential isoform expression was calculated based on the total expressed number of isoforms per gene using psiPerIsoform included in SUPPA2 together with diffSplice to compare differences in isoform expression between two conditions.
Isoform expression shifts between conditions or species were quantified using PSI values (proportion of spliced isoforms) assigned for each alternatively spliced isoform per gene as given in the equation below. We used the PSI values to calculate Shannon entropy per gene as described by Ritchie et al. (2008) 46 and used normalized values between 0 and 1 for between species comparisons as described in Kumar et al. (1986) 47
We calculated Shannon entropy values for genes expressed in control and salt treated samples for S. parvula and A. thaliana. Genes with PSI values less than 0.01 or higher than 0.99 (expected when an isoform is rarely expressed or dominates approximating zero alternative splicing for that gene) were removed from our analysis to test for isoform expression shifts. Further, gene models which were not represented by at least two isoforms were removed from the analysis.
We used a deep-neural network, SpliceAi 48 to predict genome wide splice sites for S. parvula from primary gene model sequences. The network model was trained first with A. thaliana gene models from chromosome 1 to 4 and validated with chromosome 5 gene models described in Araport11 49 . We provided 200 nucleotides upstream and downstream of a given base scanning all bases per gene in all genes models to predict whether that site is a splice site donor, acceptor or not a splice site. Model prediction was assigned a probability score between 0 and 1 for a given site with values closer to 1 representing the probability of that site being a splice site. We used a probability score of ≥0.6 for the selection of potential splice sites. We used this trained network to predict splice sites for the S. parvula genome v2. We compared the predicted splice sites to observed splice sites and identified new splice sites. If new splice sites were predicted for a gene model we had identified more than one isoform, the prediction of a novel splice site or sites for that gene was considered as one additional predicted putative isoform.
Prior to this study, the S. parvula reference gene models (v2.2) were predicted based on ab inito methods as well as RNA-seq evidence based prediction derived from non-stressed conditions 29 . To maximize the identification of transcripts that may be conditionally expressed under stress, we used 4-week-old S. parvula plants treated with multiple salts (NaCl, KCl, LiCl, and H 3 BO 3 ) that are found at high levels in its native soils 50 for PacBio Iso-Seq sequencing. We obtained 500,265 error corrected circular consensus sequences (CCS) as our primary source of sequence reads to create an isoform specific reference transcriptome and to supplement the genome-based transcript annotation for S. parvula. We identified putative full-length transcripts based on 338,812 high quality CCS reads that contained 5' and 3' primers and polyA tails (Figure 1A). Following iterative clustering, error correction, and mapping to the S. parvula reference genome, we annotated 16,828 (corresponding to 11,348 genomic loci) structurally distinct putative protein coding transcripts expressed in S. parvula tissues exposed to multiple salts (see Methods for details). This added 7,732 new protein coding transcript models to the S. parvula reference genome to provide a total of 34,582 reference protein coding transcripts (Table 1).
We were able to improve the S. parvula reference genome to include full length transcripts inclusive of 5' and 3' UTR regions with Iso-Seq reads. The average length of new Iso-Seq supported reference transcript models was greater than the corresponding length of transcript models in the S. parvula v2.2 reference genome annotation (Fig. 1B). The increase in transcript lengths was largely due to the identification of 5' and 3' UTR sequences that were previously missed in transcript model predictions in the reference genome. This refinement of reference transcript models generated UTR length distribution comparable to that of A. thaliana reference genome (Araport11) (Figure S1) and significantly increased the percentage of standard RNA-seq reads mapped to the reference transcriptome (Fig. 1C). This is expected to improve estimates of gene expression counts when using short-read RNA-seq data.
New genes previously not reported for S. parvula was added with Iso-Seq supported transcripts. The current reference S. parvula v2.2 genome includes 26,847 total protein coding primary gene models. The Iso-Seq supported transcripts mapped to 11,348 (42%) of those genes (Table 1). We additionally identified 301 (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for this preprint this version posted October 17, 2022. ; https://doi.org/10.1101/2022.10.13.512046 doi: bioRxiv preprint novel gene models that were missed (i.e. sequence present in the genome but annotation absent) in the S. parvula reference genome. For example, the putative ortholog of the Arabidopsis Magnesium/proton exchanger (MHX, similar to At2G47600) in the S. parvula genome was annotated on chromosome 4 between Sp4g29520 and Sp4g29540, using an Iso-seq based transcript model detected in this study (Figure S2). The novel transcript models further improved the reference genome annotation by adding multiple isoforms assigned to gene models, alternative transcription start and end sites for existing models, and UTR sequences (Table 1). The improved gene models, isoform specific expression, and Iso-Seq reads are available at Bioproject ID PRJNA63667.
We identified 5,911 alternatively spliced events, resulting in structurally different protein coding regions from the primary transcript models in the S. parvula genome from this study. These splice variants were categorized into intron retention, alternative 3' acceptor, alternative 5' donor, exon skipping, use of alternative first exon, use of alternative last exon, and use of mutually exclusive exons based on their frequency (Table 2). Intron retention was the most prevalent (55.2%) alternatively spliced event in S. parvula.
We observed that two or more distinctly spliced isoforms could be co-expressed in either shoots or roots (Figure S3) when multiple isoforms were checked for their expression using RT-PCR for a select set of genes.
Alternative splicing can increase the repertoire of transcripts that are available to respond to abiotic stresses more efficiently and dynamically, independent of gene copy number variation 51 . Therefore, we hypothesized that S. parvula would have a higher diversity of alternatively spliced isoforms for genes related to abiotic stress tolerance, specifically salinity tolerance, than in the less-tolerant species A. thaliana. To test this, we calculated the isoform ratio per ortholog pair in S. parvula and A. thaliana using the S. parvula reference isoforms identified in this study and A. thaliana reference isoforms obtained from AtRTD2 database 16 . To avoid missing data or lack of expression of a certain gene in mature shoots or roots in one species being inferred as lack of isoform diversity in that species, we limited our comparison to genes expressed in our study that were represented by at least one transcript model in both species. We identified 10,859 A. thaliana -S. parvula ortholog pairs that had one or more isoforms per ortholog in each species (Figure 2A; Supplementary Table 2). Among them there were 6,874 ortholog pairs showing more isoforms in A. thaliana while only 1,201 pairs had a higher isoform number in S. parvula (Fig. 2A). Ortholog pairs annotated as "Response to stress" (GO:0006950) and "Transport" (GO:0006810) had a higher isoform diversity in S. parvula, while ortholog pairs annotated under "Nitrogen metabolism" (GO:0034641) had a higher isoform diversity in A. thaliana (Fig. 2B). As a control, we examined the distribution of isoforms in all ortholog pairs and found that these distributions were not significantly different between the two species (Fig. 2B).
The genes that had a higher isoform diversity in S. parvula included some of the most highly conserved and key stress responsive genes in plants including the Na + /H + antiporter, SOS1 known for its role in excluding Na + from roots during salt stress 52 and P5CS1 that codes for delta1-pyrroline-5-carboxylate (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for this preprint this version posted October 17, 2022. ; https://doi.org/10.1101/2022.10.13.512046 doi: bioRxiv preprint 9 synthase, the rate-limiting enzyme in proline biosynthesis known for its role in oxidative and osmotic stress responses 53 . Notably, both SOS1 and P5CS1 are represented by single copy orthologs in S. parvula and A. thaliana. Five SOS1 (out of 8 detected) and 8 P5CS1 (out of 22 detected) isoforms for S. parvula were expressed at ≥ 0.5 TPM in both shoots and roots in control as well as salt treated conditions (Figure S4). The AtRTD2 database reported three SOS1 and six P5CS1 isoforms for A. thaliana 16 .
Isoform usage is less disordered in S. parvula compared to A. thaliana during salt stress Diversity and conditional expression (i.e. specificity) of isoforms can be assessed using the Shannon entropy based information theory applied to transcriptomes 54 . Stressed compared to growth optimal conditions are known to have higher transcriptome entropy and disorderdness with an increased number of alternative splice events when assessed using Shannon Entropy 46 . We hypothesized that S. parvula transcriptomes will show a smaller entropy increase in its isoform usage when transitioning from control to salt stressed treatments compared to the salt-sensitive model A. thaliana. To test if isoform usage from control to stressed conditions went through a measurable entropy transition distinctive of the species, we used RNA-seq data from root and shoot samples to quantify the isoform abundance in S. parvula and A. thaliana and calculated the Shannon entropy (see Methods). We used A. thaliana-S. parvula ortholog pairs that were represented by at least two expressed isoforms with a normalized expression ≥ 0.5 TPM per ortholog within a species to avoid incomplete comparisons due to rare isoforms difficult to quantify in one species. This resulted in a total of 1,678 and 1,592 ortholog pairs expressed in roots and shoots. Roots had 3,832 and 5,239 isoforms for S. parvula and A. thaliana while shoots had 3658 and 4431 isoforms respectively. We found that both S. parvula and A. thaliana root transcript distributions increased mean entropy in response to salt stress (Figure 3A). This is aligned with the expectation that stress conditions create higher transcript diversity, lower specificity, and more disorderdness in transcript expression compared to a stress-neutral control condition 55 . A. thaliana shoots showed a significant increase in entropy when transitioning from control to salt stressed conditions (Fig. 3A).
The change in entropy for S. parvula was less in both roots and shoots suggesting a less disordered state of isoform usage compared to the relatively stress-sensitive A. thaliana when responding to stress conditions. We observed that the isoform usage in response to salt was highly species specific. The number of ortholog pairs that showed increased or decreased isoform usage as a shared response to salt stress in both species roots (156 expressed orthologs) and shoots (142 expressed orthologs) were much fewer than those orthologs (977 in roots and 937 shoots) that had a specific usage change in one species (Fig. 3B). Orthologs that showed high isoform usage specificity (i.e. maintained or lowered entropy) in response to salt stress in S. parvula roots compared to A. thaliana were enriched in functions largely associated with salt stress (Fig. 3C). In shoots, genes that maintained isoform usage specificity under salt stress in both species were enriched in salt stress associated functions (Fig. 3C).
(which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for this preprint this version posted October 17, 2022. ; https://doi.org/10.1101/2022.10.13.512046 doi: bioRxiv preprint 10
We next examined if the differently expressed genes in response to salt stress were also subjected to changes in their isoform usage under high salinity. Supplementary Table 3 lists all genes identified as differently spliced or differently expressed. Genes that were differently expressed as well as differently spliced in response to salt stress were rare in S. parvula and A. thaliana (≤ 3%) (Figure 4A). Moreover, the shared orthologs that were differently spliced in response to salt stress between species either in roots or shoots were also low (~3%) (Fig. 4B). Multiple genes differently expressed under salt stress in A. thaliana are found to be only differently spliced in response to salt in S. parvula (Supplementary Table 3). Figure S5 further highlights the high degree of species-specific regulation in differential isoform usage in response to stress. However, there is high convergence in the enriched functions represented by differently spliced isoforms in response to salt stress in both species (Fig. 4C).
Majority of splice sites in plants are marked by GU at the 5' and AG at the 3' sites in introns 56 .
Although less common, plant genes are spliced at alternative sites termed as non-canonical splice sites and alternative splicing at non-canonical sites are associated with abiotic stress responses [56][57][58] . We investigated whether the expression of transcripts with non-canonical splice sites (Supplementary Table 4) increased under salt stress in S. parvula compared to A. thaliana. We found that S. parvula did not show any significant difference in mean expression strength between non-canonical and canonical transcripts in both roots and shoots while A. thaliana shoots showed an increased expression in transcripts that had non-canonical splice sites when treated with salt (Figure 5A).
Previous studies have reported increases in non-canonical splicing in plants under abiotic stresses 59,60 .
Therefore, we examined if usage of transcripts with non-canonical splice sites significantly increased under salt stress compared to control conditions in S. parvula differently from A. thaliana. Similar to previous reports, non-canonically spliced transcripts are less frequent than canonically spliced transcripts regardless of the condition tested (≤ 10%; Fig. 5B). However, our analysis does not find a significant increase in noncanonically spliced isoforms from control to salt treated conditions in either species (Fig. 5B).
Next, we tested if genes with non-canonical splice sites were enriched for stress associated functions in S. parvula. We found 424 genes out of 25,145 multi-exon coding genes to be enriched in non-canonical splice sites in the S. parvula genome (Fig. 5C). Some of these are notable genes associated with stress regulatory pathways (for example, PAL1, PAL2, P5CS1 and HSC70-1) (Fig. 5C). S. parvula genes enriched for non-canonical splice sites were indeed primarily enriched in stress response pathways (Fig. 5D). Further subclustering of the functional group annotated under "stress responses" (cluster C1 of Fig. 5D) showed that genes in salt/metal ion and osmotic stress were specifically contributing to this cluster.
It is likely that we may have missed to detect stress responsive isoforms expressed in S. parvula in this (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for this preprint this version posted October 17, 2022. ; https://doi.org/10.1101/2022.10.13.512046 doi: bioRxiv preprint 11 study because exhaustive searches for conditionally expressed isoforms are impractical for emerging model organisms. Therefore, we sought to employ a machine learning approach to predict alternative splicing sites in the S. parvula genome as an alternative. We applied the deep neural network, SpliceAI which is expected to yield high confidence predictions among recent tools developed to predict splice events using genomic sequences 48,61,62 . We used the known splice site information from A. thaliana chromosomes 1-4 to train the SpliceAI network and received an average precision of 0.92 when tested with A. thaliana chromosome 5 (Figure 6A). We then predicted splice sites from 26,847 S. parvula pre-mRNA sequences and obtained 214,901 splice site predictions including 114,284 novel splice sites (Fig. 6B). Twenty-six percent of splice sites previously observed were also predicted using SpliceAI and we found 7,302 genes with at least one newly, predicted isoform. Prediction probability scores were highest for splice sites within the gene compared to those in the first and the last introns (Figure S6).
With the current analysis, we have identified 16,061 potential protein coding isoforms (observed or predicted) for 9,033 genes in the S. parvula genome (Supplementary Table 5). Interestingly, stress and transport associated functions are enriched among those genes that are observed or predicted to have more than one isoform (Fig. 6C). Stress and transport related functions deduced from GO annotations account for 35% of genes that are alternatively spliced in the S. parvula genome.
Alternative splicing allows genes to acquire new functions independent from gene duplications and promoter evolution. Previous studies have shown that duplicated genes are enriched in stress associated functions in S. parvula and other extremophytes facilitating their stress adapted lifestyles more than in stresssensitive sister species 29,63,64 . However, extremophyte gene diversity represented by alternatively spliced isoforms is underexplored 65 . Certain genes are regulated only at the alternative splicing level with no change at the gene expression level that have led to the increasing recognition of the importance of isoform specific reference transcript datasets in gene expression studies 66 . In this study we examined the possibility of diversifying gene functions through alternative splicing and specially focused on isoforms differently used during salt stress in one of the leading model extremophytes 26 .
Schrenkiella parvula and A. thaliana genomes have similar gene numbers (~27,000) and similar genome sizes (~120 MB) 29 . A recent study that explored alternative splicing in A. thaliana using full length transcript sequencing based on Iso-Seq reports the discovery of isoforms in similar proportions to our study with intron retention being the most common alternative splicing event 67 (Table 2). This suggests that S. parvula is not an exception in highly increased or decreased transcript diversity through alternative splicing although the recorded number of isoforms for the model plant through aggregate studies using multiple tissues, developmental stages, and treatments are much higher (Zhang et al., 2017). Given the genomic similarities (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.
between S. parvula and A. thaliana, their transcriptome adjustments with differential splicing in response to salt stress were remarkably distinct from one another when an identical salt treatment was given to mature plants (same age and tissues tested in both species) (Figs. 4 andS5).
In support of our hypothesis that extremophytes would diversify their response to stress via alternative splicing in selected gene groups, we observed that S. parvula orthologs had a higher number of isoforms compared to A. thaliana in genes associated with stress and transport functions (Fig. 2). Stress and transport functions were also enriched among duplicated genes in S. parvula compared to A. thaliana 68 . We found that differently expressed genes and genes that showed differential isoform usage were largely mutually exclusive within species as well as in one-to-one ortholog pairs between S. parvula and A. thaliana (Figs. 4 andS5).
Further, when we combine both observed and predicted splice sites in the S. parvula genome, the potential protein coding isoform pool is enriched in functions associated with stress tolerance (Fig. 6). These observations together indicate that genes expressed in response to stress are highly diversified and nonoverlapping in their mode of function, but converge on common functions associated with stress tolerance in S.
parvula. Therefore, our study provides a novel resource for assessing functional significance of stress tolerance genes in the extremophyte model. It allows selection of target genes that could be tested at the isoform level when expression modulation via promoter modifications or single gene-knockouts of essential genes do not offer optimal methods to test novel gene functions contributing to stress tolerance.
Our current study in agreement with a previous study on A. thaliana have shown that most differently spliced genes were not differently expressed in response to salt stress representing an independent layer of gene regulation in response to stress 21 . Compared to animals, plants tend to use alternative splicing biased to environmental stress responses more than for tissue-specific responses 69 . Multiple studies have reported specific associations of alternative splicing and environmental stress in plants 22,57,59,70,71 . However, fewer studies have examined the presence of non-specific alternative splicing leading to increased number of differently spliced isoforms under abiotic stress 12,72 . Additionally, components of the spliceosome are differently expressed leading to differential splicing of target genes in A. thaliana during stress conditions 60 .
In animals, stressed conditions are reported to have increased amount of alternatively spliced isoforms with high non-specific expression, thus creating a higher level of disorderdness in isoform expression 46 which can be quantified using Shannon entropy 73,74 . We predicted that plants will show a similar trend in increased disorderdness in isoform expression at the transcriptome level during stress conditions. Furthermore, we expected to see a smaller change in entropy in the extremophyte when transitioning to a salt treated condition compared to the stress-sensitive species. Indeed, this prediction was supported by the shoot transcriptomic response we observed for S. parvula and A. thaliana (Fig. 3). Notably, the genes that shifted to lower entropy values representing shifts to specific isoform in their expression specificity under stress were enriched for stress associated functions in both roots and shoots in S. parvula (Fig. 3). Our study cannot test if the tendency (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for this preprint this version posted October 17, 2022. ; https://doi.org/10.1101/2022.10.13.512046 doi: bioRxiv preprint 13 to increase transcriptome disorderdness via less specifically expressed isoforms per gene is indicative of aberrant splicing under stress. Yet, the comparison between S. parvula and A. thaliana suggests that the extremophyte is more prepared to respond to salt stress by specific isoforms mostly expressed for stress associated genes.
In conclusion, this study provides a novel resource for a leading extremophyte model and expands our knowledge on the ability to respond to stress via differential isoform usage independently from differential gene expression. Stress associated functions were enriched among genes observed or predicted to have multiple isoforms in S. parvula; one-to-one orthologs where S. parvula has a higher number of isoforms than A. thaliana; genes that showed differential isoform usage in response to stress in S. parvula; S. parvula genes that were enriched in non-canonical splice sites; and S. parvula genes that maintained or lowered their disorderdness by expression of specific isoforms under stress. These findings contribute to how we understand stress tolerance evolved in an extremophyte. Differential isoform usage offers a complementary path to increase the coding potential of the S. parvula genome that cannot be fully explained by gene duplication or promoter evolution alone. Future studies on other extremophytes exploring isoform diversity will facilitate the identification of convergent traits in isoform usage evolved in stress-adapted plants. Such a resource will be influential in deducing diverse stress responsive networks and identifying transferable stress responsive genes into crops. Table 1. Summary of S. parvula V2 genome updated with transcript models supported by Iso-Seq
Genes supported by an Iso-Seq based transcript model 11,348
Total transcript models identified from Iso-Seq 16,828
Gene models identified with at least one new isoform 7,028
New genes identified from Iso-Seq 301
Gene models supplemented with UTR information using Iso-Seq 11,348
Genes with alternative splicing identified using Iso-Seq 3,470
Genes with alternative starts identified using Iso-Seq 4,760
Genes with alternative ends identified using Iso-Seq 4,756
Total number of protein coding transcripts annotated in the genome 34,582 [C] Functionally enriched processes represented by ortholog pairs in distinct categories of entropy shifts. A node in each cluster represents a gene ontology (GO) term; size of a node represents the number of genes included in that GO term; the clusters that represent similar functions share the same color and are given a representative cluster name and ID; and the edges between nodes show shared genes between functions. All clusters included in the network have adj p-values ≤0.05 with false discovery rate correction applied.
(which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for this preprint this version posted October [A] Expression distribution of transcripts that contain only canonical splice sites and transcripts with at least one non-canonical splice site in roots and shoots of S. parvula (left panel) and A. thaliana (right panel). Asterisks indicate significant difference of expression distributions between control and salt treated condition measured by two-sided t-test at p-value ≤ 0.05. [B] Number of expressed non-canonically spliced transcripts as a % out of total transcripts expressed in S. parvula and A. thaliana in response to salt. Significant differences between control and stress conditions were tested using Fisher's exact test. [C] S. parvula genes that were enriched in non-canonical splice sites. The y axis shows the -log 10 p-value for a test of excess of non-canonical splice sites computed using a binomial test, where the probability of enrichment is calculated as the total non-canonical splice sites divided by the total number of splice sites per gene, ordered in the chromosomal order (x-axis) for the S. parvula genome. Genes with a high enrichment for non-canonical splicing are labeled. Red line indicates thelog 10 p corresponding to adjusted p-value of 0.05. [D] Functional processes enriched in genes detected to be noncanonically spliced in [C]. A node in each cluster represents a gene ontology (GO) term; size of a node represents the number of genes included in that GO term; the clusters that represent similar functions share the same color and are given a representative cluster name and ID; and the edges between nodes show the connectivity of genes between functions. All clusters included in the network have adj p-values ≤0.05 with false discovery rate correction applied. More significant values are represented by darker node colors. The right panel shows the sub-clustered functions represented by the largest cluster C1 in the left panel.
(which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for this preprint this version posted October 17, 2022. ; https://doi.org/10.1101/2022.10.13.512046 doi: bioRxiv preprint [C] Functional processes enriched in genes observed and predicted to have more than one isoform in the S. parvula genome. A node in each cluster represents a gene ontology (GO) term; size of a node represents the number of genes included in that GO term; the clusters that represent similar functions share the same color and are given a representative cluster name; and the edges between nodes show the connectivity of genes between functions. All clusters included in the network have adj p-values ≤0.05 with false discovery rate correction applied.
(which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.The copyright holder for this preprint this version posted October 17, 2022. ; https://doi.org/10.1101/2022.10.13.512046 doi: bioRxiv preprint