Polypedilum vanderplanki or the sleeping chironomid is a dipteran in the family Chironomidae (non-biting midges). Northern Nigeria and Uganda are the main habitats of the species. The first time it was discovered by Dr. F. L. Vanderplank, who noticed that larvae of the species can survive in dried mud. The main characteristics of chironomids and effects of high temperature on larvae first time were described by H. E. Hinton in 1951. [19]

Figure 2.1. Schematic illustration of Polypedilum vanderplanki life cycle, included 5 stages of anhydrobiosis used in the research. Figure (with modification) was taken from [21, 51].

Anhydrobiosis is the most studied form of cryptobiosis Cryptobiosis - a physiological state in which metabolic activity is reduced to an undetectable level without disappearing altogether entered in response to adverse environmental conditions. and occurs in situations of extreme desiccation. (from Greek `life without water') It was observed in many microorganisms, plant seeds, some plants, invertebrate animals and at least one insect -- Polypedilum vanderplanki.[10, 52, 56]

However, the sleeping chironomid exhibits desiccation tolerance only on the larvae stage. In anhydrobiotic condition larvae can live more than 17 years and be exposed to temporary high (103є) or low (-190є) temperatures, radiation (7000 Gy) or immersing into 100% ethanol. [16, 57]

The process of desiccation takes around 48 hours, while rehydration starts in 30 minutes after contact with liquid water and continues up to 24 hours. (Figure 2.1) During these short periods, the organism has to go through complicated biochemical processes determining further survival of larva. [7, 19]

Trehalose, hemolymphatic disaccharide of insects, plays the key role in successful desiccation. It is the well-known fact that trehalose works as a compatible solute in an organism with desiccation tolerance, such as nematodes, tardigrades and Artemia cysts. Due to chemico-physical characteristics, this disaccharide can be accumulated in a high quantity without harm to organisms, replacing water and protecting biomolecules from degradation. [13] There are 3 theories describing the trehalose protection mechanism. The first theory suggests that trehalose forms hydrogen bonds inside and outside the bilipid membrane to maintain its integrity instead of water. The second hypothesis is based on the fact that dehydration causes the formation of `trehalose glass', which fills the space originally occupied by water, and limits the physical and chemical interactions between biomolecules and their mobility, thereby preventing their damage. The third is based on the assumption that trehalose retains a thin layer of water around the proteins, which prevents their denaturation during drying [8, 9, 48]. In the sleeping chironomid trehalose synthesis grows rapidly when the water content in the organism decreases to 75%. Presumably, the increased osmotic pressure and changes of ionic concentration are the factors inducing the trehalose synthesis. [13, 58]

In addition to trehalose, Late Embryo Abundant (LEA) proteins are also involved in the protection of cells from damage during desiccation. Initially, they have been found in plant seed and shown its accumulation in response to dehydration and changes in ion concentration. The protein changes its tertiary structure from random to coiled-helix structure during desiccation and, hypothetically, binds with intracellular ions or with each other for cell protection.[24] gene splicing cell polypedilum

Moreover, Heat shock proteins (HSPs) also work as protein protectors during desiccation in P. vanderplanki. HSP belongs to the chaperones family and produced in the cell under stress environments such as extreme temperature or UV light.[54] In the sleeping chironomid several HSPs increased their expression during dehydration were found (Pv-hsp90, Pv-hsp70, Pv-hsc70, Pv-hsp60). Thus, during the dehydration-rehydration cycle in P. vanderplanki some heat shock proteins protect proteins from denaturation and refold denatured proteins. [15]

Along with trehalose, LEA and HSPs, antioxidants, small redox proteins, protein-repair methyltransferases, hemoglobins, and aquaporins also participate in the successful course of anhydrobiosis. [16]

Polypedilum vanderplanki continues to be one of the favorite model organisms to study the anhydrobiosis process. Scientists keep exploring the genome, transcriptome, and proteome of the sleeping chironomid, looking for changes between stages of the dehydration-rehydration cycle and regulatory pathways brought to them [15, 18, 39].

1.2 Mechanism of splicing

Since 1977 science world considers that eukaryotic genes contain a lot of interruptions and mRNA is not just a copy of DNA as an indisputable fact. They consist of protein-coding regions (exons) and non-coding regions (introns). Around 15 years later Phillip Allen Sharp and Richard J. Roberts discovered the RNA splicing, the process of introns removing and transformation of pre-mRNA into mRNA.[1, 6, 49] However, not only protein-coding genes go through RNA processing steps. The majority of long non-coding RNA also undergo splicing. [29]

The key roles in the RNA splicing process play splicing sites - short conserved sequences at intron ends. Spliceosome, a complex molecular machine consisting of multiple RNAs and protein, recognizes and bind splicing sites on newly synthesized RNA. There are three splice sites in an intron. First, donor 5' splice site (5'SS), it locates close to the 5' end of the intron and contains conservative dinucleotide - GU. Second, acceptor 3' splice site (3'SS) places near the 3' end of the intron and is represented by almost invariant dinucleotide - AG. Upstream of 3' SS is a region with high C and U content, polypyrimidine tract. Third, around 50 nucleotides upstream from the acceptor site lie the branchpoint (BP), which included an adenine nucleotide. (Figure 2.2)[50]

Figure 2.2. Logos of splice sites (a) 5' splice site (b) branch point site (c) 3' splice site. [30]

From the chemical side of the question, splicing is two sequential phosphoryl-transfer reactions (transesterification). During first, branching, 2' hydroxyl of adenosine in the branchpoint attacks a phosphate at 5' SS, thereby splitting the upstream exon and 5' end of the intron and forming an intermediate stage with `lariat'. In the second reaction, exon ligation, the new formed 3' hydroxyl group on 5' end attacks phosphate at 3' SS to merge the exons and remove the intron. (Figure 2.3) However, easy on the first glance the process becomes more complicated if taking into account the spliceosomal complex.[20]

Figure 2.3. The schematic view of the splicing process. (top) Two sequential phosphoryl-transfer reactions[14] (bottom) Spliceosome assembly and transformation complex, where U1,U2,U4,U5 and U6 are the spliceosomal small nuclear ribonucleoproteins. [30] Abbreviation: SS, splice site.

A spliceosome is a ribonucleoprotein enzyme, which catalyzed the splicing reaction. It is composed of snRNAs and over 70 associated proteins. Major spliceosome contains 5 spliceosome subunits: U1, U2, U4, U5, and U6. They in combination with proteins form complexes in different conformations and generate an active enzyme site. (Figure 2.3)

The first spliceosome complex ('A') forms when U1 and U2 bind with 5' splice site and branch point, respectively. Three subunits, U4, U5, and U6, generate triple snRNP, bind with the complex A and transform it into the precatalytic spliceosome ('B'). However, complex B is not stable and active, thus ATPases Prp28 and Brr2 dissociate U1 and U4 snRNA, respectively. A new B^act complex contains U2, U6, and U5 catalyze the first phosphoryl-transfer reaction to form the C complex. The ATPase Prp16 remodels spliceosome again and prepares it for the `intron ligation' reaction. The U5 subunit aligns flanking exons and finally, the intron is ligated, flanking exons merge and mRNA releases. [14, 59]

There are two potential mechanisms of splice site recognition called intron and exon definition. (Figure 2.4) The first is based on intron identification and places the splicing basal machinery across introns. The length of the intron, in this case, should be quite small (200-250 nt). This way could be more ancient and spread in organisms with short introns, for example, yeast. In higher metazoa, particularly vertebrates, the opposite situation happens. There are short exons (~300 nt) and long introns. Thus exon definition is more logical and the basal splicing machinery is placed across exons. However, the later spliceosome actions are aimed at the intron region.[45, 60]

Figure 2.4. Schematic illustration of intron and exon definition processes, where U1 and U2 are the spliceosomal small nuclear ribonucleoproteins and ss - splice site.[5]

Alternative splicing depends on splice sites' selection by spliceosome. There are two models explaining how the transcriptional environment can affect this choice.

The `recruitment model' suggests that alternative splicing is regulated by the direct or indirect recruitment of trans-elements to pre-mRNA. The `kinetic model' proposes that the rate of RNAPII elongation affects the process of splice sites' selection. [12]

The RNA polymerase II speed. According to the kinetic model, the changes in the rate of RNAPII complexes can affect alternative splicing regulation. Under slow transcription spliceosomes have enough time to assemble on weak splice sites that results in increased inclusion of cassette exons. Accordingly, faster elongation results in a greater number of splice sites available simultaneously, that allow spliceosome to select stronger one.

Differential splicing analysis. Differential splicing analysis was carried out by SAJR.[38] It splits the genome into segments, regions between 2 splicing sites, count the reads which overlapped or not the segment and inclusion ratio (psi, ??). (Formula 1 and Formula 2 for alternative first and last exons )

Also, it identifies `Constitutive' and `Alternative' segments, existing in all transcripts or only in part of them, respectively. (Figure 3.1)

, (1)

, (2)

Figure 3.1. A schematically step-by-step explanation of test and control segments' definition.

Where i - number of inclusion reads, e - number of exclusion reads, lr - length of a read, ls - length of a segment.

When all measurement was finished, segments with significantly changed inclusion ratio between stages of the dehydration-rehydration cycle were found by building and fitting GLM model in R SAJR package and Benjamini-Hochberg (`BH') correction. The threshold for the p-value was set up on 0.05. The segments were divided on up- and down-regulated by the delta PSI between stages with a threshold in 0.1. In further operations, the up- and down-regulated segments will be taken as test dataset, while resting of alternative and all constitutive segments as control datasets.

Differential expression analysis. Differential expression (DE) analysis was carried out by edgeR package in R based on genes count from SAJR.[40] Data were filtered by the expression, normalized and differentially expressed genes were figured out by GLM model and `Benjamini-Hochberg' correction. The threshold for the p-value was set up on 0.05. (Table 6.2)

Functional annotation. The most probable cDNA and corresponding proteins were extracted from the genome by `Transdecoder' and uploaded to `Interproscan' for further functional annotation of proteins. [17]The program scan sequences and find matches in the following databases: TIGRFAM, SFLD, ProDom, Hamap, SMART, CDD, ProSiteProfiles, ProSitePatterns, SUPERFAMILY, PRINTS, PANTHER, Gene3D, PIRSF, Pfam, Coils, MobiDBLite. [42]

GO enrichment analysis. So far as a functional annotation with GO id has been obtained, GO enrichment analysis was done. Genes were filtered by several parameters.

Figure 3.2. All used combinations of test dataset and background for GO enrichment analysis.

They must have GO id, count per million more than 1 on all stages of the dehydration-rehydration cycle, and contain at least one internal exon. The genes were separated on geneset by the GO id and used on the background. On the figure 3.2 below showing all variants of comparisons. Significantly enriched gene sets were identified by Fisher's exact test and `BH' correction with a threshold of 0.05.

Strength of the splicing sites. Position weight matrix (PWM) was built for splice sites and polypyrimidine tracts based on the constitutive segments by R package `Biostrings'. For both test and control datasets scores were calculated and compared by Student's t-test with the threshold for p-value in 0.05. (Figure 3.3)

Figure 3.3. Illustration of introns' selected zones for exploring the strength of the splicing sites.

One more mark of polypyrimidine tract strength is polypyrimidine enrichment of the regions. Then bigger is the purine content in the polypyrimidine tract than more reducing the affinity with U2AF. [22] CT content was found by the simple formula.

,(3)

Where C - number of cytosines, T - number of thymines, sl - sequence length.

Motif enrichment analysis. As far as motif enrichment analysis of complete introns is a long process, in the research it was done only on certain intron regions. Based on the intron length histogram, the minimum length was set on 100 nucleotides. Then without splicing sites and polypyrimidine tract, it will shorten up to 75 nucleotides. Thus test datasets contained cassette exons flanked introns not shorter than 100 nt. In the controls were selected the segments from the same genes surrounded not shorter introns.

These sets were pairwise compared to the existing of the Drosophila motifs from the CISBP-RNA database. [46] The motif enrichment analysis carried out with the `AME' tool from web-server MEME.[41] It was run with the following settings:

ame --verbose 1 --oc . --scoring totalhits --method fisher --hit-lo-fraction 0.25 --evalue-report-threshold 10.0 --control exn_upstream.fasta d_se_upstream.fasta db/CISBP-RNA/Drosophila_melanogaster.dna_encoded.meme

It uses Fisher's exact test and Bonferroni correction with threshold 0.05 for finding the significantly overrepresented motifs.

Orthologs' search. All known RNA binding proteins from the CISBP-RNA database for found enriched motifs around test segments were downloaded and concatenated with Polypedilum vanderplanki proteome from Transdecoder. This collection of amino acid sequences was uploaded as a BLAST database and ran in protein BLAST in all-against-all mode. For each motif were selected the best hit one or several proteins by E-value.[3]

3. Results and Discussion

3.1 Differential splicing analysis

The RNA-seq data consisted of five stages of dehydration/rehydration cycle (fully hydrated, 24 hours of dehydration, 48 hours of dehydration, 3 hours of rehydration and 24 hours of rehydration, two biological replicates for each stage) were mapped to the reference genome. The reads were counted, and differential splicing analysis was performed using SAJR pipeline [38]. Multidimensional scaling (MDS) plot confirms data consistency. As the PvD48 samples are most remote from PvD0 and PvR24 are the closest to it. (Figure 4.1)

Figure 4.1.Inclusion ratio MDS plot

The 1322 segments were identified as significantly changing their inclusion between fully hydrated stage and at least on stage of dehydration/rehydration cycle. (GLM, quasi likelihood ratio test, p-value < 0.05, |dPSI| > 0.1). 688 and 528 segments were up- or down-regulated at one or more stages of rehydration/dehydration cycle compared to fully hydrated form. For 106 segments change direction was different at different stages. (Table 4.1, Figure 6.1).

Table 4.1. The number of up- and down-spliced segments on different stages of the dehydration-rehydration cycle.

For down-regulated segments PSI distribution is almost uniform before desiccation and strongly skewed to zero on the vitrification stage (d48), while the up-regulated segments exhibit the opposite behavior (Figure 4.2).

However, the majority of segments demonstrated moderate amplitude of PSI change under dehydration (Figure 4.2.C).

The delta PSI mostly lay within 0.1 and 0.2, which means that the inclusion of the region reduced or increased on average 10-20% from the normal condition.

The inclusion of alternative segments' various types changes in different ways between stages of the dehydration-rehydration cycle.

Cassette exons are more often excluded from the transcript during desiccation with the peak of dropout on 48 hours, while intron retention, oppositely, increases. Alternative donor and acceptor segments also begin to be included in transcripts during dehydration and even more at the beginning of rehydration. (Table 4.1, Figure 4.3)

Thus, cassette exons' PSI decrease during desiccation, while intron retention, alternative donor and acceptor segments start to include in transcripts.

Figure 4.2. (A and B) The distribution of down- or up-regulated segments' inclusion ratio in normal active state (d0) and 48 hours of dehydration (d48). (C) The distribution of delta PSI (dPSI) for down-regulated (negative_dPSI) and up-regulated (positive_dPSI) segments at 48 hours of dehydration.

Figure 4.3. The distribution of different types of segments between anhydrobiosis stages. Abbreviations: d24 - 24 hours of dehydration; d48 - 48 hours of dehydration; r3 - 3 hours of rehydration; r24 - 24 hours of rehydration; pos_dPSI - up-regulated segments; neg_dPSI - down-regulated segments.

Using hierarchical cluster analysis with correlation distance we clustered desiccation-related segments into 5 clusters (Figure 4.4). It separated the segments into groups by patterns of their inclusion ratio changes during rehydration/dehydration cycle. In three of them segments reduce their inclusion during anhydrobiosis with a minimum PSI at 48 hours of dehydration (clusters 1 and 3) or 3 hours of rehydration (cluster 5). Segments from remaining clusters demonstrate increase of their inclusion in comparison with active larva state (clusters 2 and 4).

Figure 4.4. Cluster analysis of statistically significant differentially spliced segments.

Subsequently, all analyzes were performed exclusively on cassette segments. Firstly, because intron retention (IR) is traditionally overlooked as transcriptional noise and the methods for IR accurate detection is limited and just starting to develop. [32] Secondly, there are more cassette exons than alternative acceptors or donors in the dataset. In addition, cassette exons are predominant type AS and bringing the greatest functional changes for mRNA.[11]

3.2 Splicing sites' strength

One of the widespread observations about splicing is that alternative splicing sites are weaker than constant [28]. It causes the necessary outside intervention on alternative splice sites' selection. The strength of the splicing sites in the control and test sets was measured based on PWM and compared with each other. Both up- and down-regulated cassette exons have significantly weaker splice sites than constitutive exons from the control set (T-test, p-value < 0.05). (Figure 4.5A-B)

Figure 4.5. Pairwise comparisons of controls and test segments' (A) acceptor splice site (B) donor splice site (C) polypyrimidine tract. (D) Logo plots for acceptor, donor splice site and polypyrimidine tract. (E) CT content in polypyrimidine tract. Abbreviations: exn - constitutive segments; alt - alternative segments, but not differentially spliced; pos_dPSI - up-regulated segments; neg_dPSI - down-regulated segments. The statistically significant comparisons shown by red arrows.

However, all analysed sites are from major spliceosome GT-AG introns (Figure 4.5D). In addition, the polypyrimidine tract upstream of alternative segments is less pyrimidine enriched. (Figure 4.5C,E) The purine content in the PPT reduces affinity with U2AF and, consequently, contributes to skipping a segment.[22]

Summarizing, in P. vanderplanki all common observations for other species splicing are kept, and should exist the outside factors underlie the observed changes in alternative splicing during anhydrobiosis.

3.3 Intron length distribution

For further motifs analysis, the knowledge of intron length distribution is necessary, because it gives an understanding of the minimum threshold for flanking introns length. Distribution of intron length exhibits clear peak around 60 nt, majority (66%) of P. vanderplanki introns are short with length from 50 to 100 nt (Figure 4.6). Similar distributions were found in several other model organisms with relatively short genomes, such as Drosophila melanogaster and Caenorhabditis elegans.(Figure 4.6)[33]

Figure 4.6. Distribution of introns' length (left) in P. vanderplanki (right) in S. cerevisiae (S.ce.), C. elegans (C.el.), D. melanogaster (D.me.), A. thaliana (A.th.) and Homo sapiens (H.sa.).[33]

All cassette exons have statistically significantly longer flanking introns than constitutive segments (Wilcoxon test, p-value < 0.05).

This suggests the abundance of the exon definition mechanism in cassette exons. Consequently, the regulation elements are most probably located in introns (Figure 4.7, Table 6.3).

Figure 4.7. The distribution of introns' length splitted by datasets, where exn - constitutive segments; alt - alternative segments, but not differentially spliced; positive_dPSI - up-regulated segments; negative_dPSI - down-regulated segments.

Down-regulated segments have even longer flanking introns than up-regulated or alternative but not desiccation-related cassette exons (Wilcoxon test, p-value < 0.05).

From figure 4.7 it became obvious that included segments are closer to other alternative ones, while excluded cassette exons are much more frequently surrounded by long introns.

The following analysis was focused on cassette exons surrounded by introns longer than 100 nt.

3.4 Motif analysis

Motif enrichment analysis. Up- and down-regulated segments flanked by introns longer than 100 nt on both sides were used as test sets in motif enrichment analysis.

We formed two control sets of segments - constitutive and alternative but not differentially spliced exons. They were selected from the same genes as tests and with long introns (> 100 bp) as well.

Motif enrichment analysis resulted in the 10 motifs enriched in the test dataset compared to one of controls (Table 6.4).

Figure 4.8. Log Fold Change during anhydrobiosis of differentially expressed genes coding the RBPs binding the overrepresented motifs.

Nine of them overrepresented in down-regulated sequences and one in up-regulated. There are 71 D. melanogaster RNA-binding proteins associated with at least one of these 10 motifs. Using BLAST we found 11 P. vanderplanki closest homologs for 58 of them.

Nine of them were differentially expressed during rehydration/rehydration cycle, in particular MSTRG.614, MSTRG.4336, MSTRG.12277 and MSTRG.15360 changed their expression more than twice during anhydrobiosis. (Figure 4.8)

Domain architecture and self-regulated segments. All proteins encoded by found homologs contain RNA recognition motifs (RRM). However, the one coded by MSTRG.614 also has the ZZ domain of cytoplasmic polyadenylation element-binding protein 1. (Figure 4.9)

Figure 4.9. Motifs logos and domain architecture of corresponding P. vanderplanki RBPs. Abbreviations: RRM - RNA recognition motif; 2m13 - the zz domain of cytoplasmic polyadenylation element binding protein 1.

Moreover, the genes MSTRG.4336 and MSTRG.10095 have alternatively spliced regions surrounded by their own recognition motifs. Thus these genes could be involved in self-regulation (Figure 4.10).

The self-regulated segment in MSTRG.4336 is a 15 nucleotide sequence located in the middle of the RRM region.

It doesn't affect a lot on amino acid sequence, however, its inclusion ratio reduces on 30 % between active larva condition and 48 hours of dehydration. At the same time, the gene expression decreased (Figure 4.10).

So this region could participate in mechanisms of preventing gene silencing or repressing the splicing factor activity, as far as it locates in the middle of RRM and the hosting gene expression decreasing.

Figure 4.10. Domain architecture, inclusion ratio and gene expression changes in genes with self-spliced segments - (A) MSTRG.4336 and (B) MSTRG10095.

The segment in MSTRG.10095 behaves in another way. It is a 140 bp sequence at the beginning of the protein and causes the stop codon if excluded. Its inclusion ratio has not changed dramatically during desiccation, however, it has increased rapidly during rehydration.

The same pattern is traced in expression, it decreased in dehydration and increased in rehydration. (Figure 4.10) Probably the splicing out of the segment led to the generation of premature termination codons (PTCs) and transcript degradation by nonsense-mediated mRNA decay (NMD). (Figure 4.11)

Figure 4.11. Schematic illustration of possible ways of splicing regulation by MSTRG.10095 encoded RNA-binding protein (M10095).

Involvement in molecular functions and biological processes. One of AME outputs is the list of true positive (TP) sequences, the segments from the test dataset with the overrepresentation of the motif in flanking introns. We calculated motif similarity based on the overlap of introns with motif matches and used this measure to split the motifs into three groups by their positive segments' intersections (Figure 4.12).

The Tra2 motif is remote from others because it is the sole motif overrepresented in the up-regulated cassette exons during dehydration. The Elav, Fne and Rbp9 also form a separate group that corresponds to motifs' nucleotide composition similarity. (Figure 4.12)

Figure 4.12. (left) Distance between groups of segments flanked by introns enriched by the corresponding motif. (right) Intersection the clusters of genes containing these segments.

The rest 5 motifs generate the third cluster, however, the Orb2 RRM stands out. The similar pattern keeps in the intersection of genes hosting the positive sequences (Figure 4.12).

Figure 4.13. (left) Molecular functions and biological processes which involve the genes with positive segments. (right) The results of GO enrichment analysis.

Figure 4.13 demonstrates molecular functions or biological processes enriched in genes containing exons that have a match of one of these 9 motifs in adjacent introns. Thus the proteins encoded by nine genes of P. vanderplanki (MSTRG.1051, MSTRG.614, MSTRG.4336, MSTRG.6810, MSTRG.7874, MSTRG.10095, MSTRG.12277, MSTRG.9458, MSTRG.15360) regulate the segments exclusion from genes that participate in such processes as mRNA processing, protein binding, signal transduction, phosphorylation, cytoskeletal protein binding and actin filament binding.

All these processes participate not only in alternative splicing auto- and cross-regulation but also in desiccation resistance in P. vanderplanki. [2, 25, 39, 44]

Conclusion