Integrative analysis of gene expression and alternative splicing in microalgae grown under heterotrophic condition

Abstract

Heterotrophic cultures are the most effective approach to overcome low growth rate challenge in the most commercial microalgae. However, the mechanism through which heterotrophic condition regulates algae metabolism are not completely clear. Alternative Splicing (AS) is a common posttranscriptional process by which transcriptome and proteome plasticity increases at different environmental conditions. To identify and characterize of AS events in Auxenochlorella protothecoides microalga grown in autotrophic and heterotrophic, RNA-Seq data were analysed. We found that AS increased with the transition from autotrophic to heterotrophic condition. 705 and 660 differentially expressed (DEG) and spliced (DAS) genes were identified for A.protothecoides was transferred from autotrophic to heterotrophic condition, respectively. Moreover, there was slight coverage between DEG and DAS genes. Furthermore, functional analysis showed that the DAS genes are most frequently related to ion binding and stimulus response. The results also indicated that prevalence of Intron retention is associated with down-regulation of the genes involved in carotenoid biosynthesis. This study provides valuable insights into transcriptional and posttranscriptional plasticity of microalgae during growth mode change.

About the Authors:

Bahman Panahi

Roles Conceptualization, Data curation, Formal analysis, Validation, Writing – original draft, Writing – review & editing

* E-mail: [email protected], [email protected]

Affiliation: Department of Genomics, Branch for Northwest & West region, Agricultural Biotechnology Research Institute of Iran (ABRII), Agricultural Research, Education and Extension Organization (AREEO), Tabriz, Iran

ORCID logo http://orcid.org/0000-0001-8523-994X

Mohammad Amin Hejazi

Roles Writing – original draft, Writing – review & editing

Affiliation: Department of Food Biotechnology, Branch for Northwest & West region, Agricultural Biotechnology Research Institute of Iran (ABRII), Agricultural Research, Education and Extension Organization (AREEO), Tabriz, Iran

Introduction

Heterotrophic cultures are the best alternatives to overcome low growth rates challenge in the most commercial microalgae [1]. In heterotrophic growth mode, organic substrate as sole carbon and nitrogen source are used [2]. Low operating cost and simplicity are the primary advantages of the heterotrophic cultivation. Moreover, heterotrophic cultivation provides an effective approach to produce different metabolites in industrial scale [3].

Although the biological responses of microalgae in heterotrophic condition have been surveyed, the underlying molecular mechanisms are not well known. It has been reported that heterotrophic growth mode affects carbon and nitrogen metabolisms [4]. Moreover, photosynthesis, carotenoids production, respiration, as well as protein synthesis are affected by changing of growth mode [5, 6].

Alternative Splicing (AS) serves a transcendental role in increasing coding capacity [7]. The functional impacts of AS fall into two main categories: protein-level alterations and transcript-level modifications [8]. It has been found that AS can act as a regulatory mechanism during developmental processes or in response to environmental conditions [9]. In addition, changing environmental condition has been shown to influence the splicing pattern [10]. Plants even seem to regulate their transcriptome post-transcriptionally in response to quickly changing environmental conditions by using alternative splicing mechanisms [11].

Like in higher plants, AS has a fundamental role in increasing biological complexity in algae. Prior study based on expression sequence tags (EST) analysis indicates that about 2.9% of all genes in Volvox carteri undergo alternative splicing [12], which is much lower than higher plants [8]. Identification of AS events based on ESTs analysis is mainly depending to the ESTs number, it is likely that with the advent of high throughput next generation sequencing (NGS) technology and development of new analytical pipelines more details about functional impacts of AS events will be deciphered [13].

Despite NGS technology has been widely applied to dissect of AS change during environmental and nutritional stresses in plant species [9; 14], it remains largely unknown whether or how AS is involved in microalgae responses to different growth condition.

To investigate the changes of AS during transition from autotrophic to heterotrophic growth condition, a set of throughput RNA sequencing data of Auxenochlorella protothecoides microalga growing at autotrophic and heterotrophic condition was retrieved from the European Nucleotide Archive (ENA) database. Then, differentially expressed (DEGs) and spliced (DAS) genes between autotrophic and heterotrophic conditions (at five time points, including hour 3, hour 6, hour 12, hour 24, and hour 72) were compared. DAS genes were mostly enriched in cellular metabolic process, nitrogen compound metabolic process, and biosynthetic process. Results of current study also demonstrated that independent mechanism of gene regulation implies by AS during transition to the heterotrophic growth condition.

Methods and material

Transcriptome data of Auxenochlorella protothecoides

RNA-seq data of Auxenochlorella protothecoides which were collected at heterotrophic growth mode (accession number: PRJNA484804) [4], were retrieved from the European Nucleotide Archive (ENA). The data contained raw paired ends data of cells grown at autotrophic growth condition (AG) and heterotrophic growth condition (HG). In AG, A. protothecoides cells were cultured in a medium containing 0.7 g KH2PO4, 0.3 g K2HPO4, 0.3 g MgSO4.7H2O, 0.3 mg FeSO4.7H2O, 0.01 mg vitamin B1, and 1 mL A5 trace mineral solution adding 5 g/L glycine under illumination of 2000 lux. The cells were grown in 100 mL flasks on a shaker set at 220 rpm and 28 ±1˚C [4]. In HG, cells were grown in a medium containing 0.7 g KH2PO4, 0.3 g K2HPO4, 0.3 g MgSO4.7H2O, 0.3 mg FeSO4.7H2O, 0.01 mg vitamin B1, and 1 mL A5 trace mineral solution adding 30 g/L glucose and 0.5 g/L glycine under dark condition.

Two biological and technical replicates were used for both autotrophic and heterotrophic condition and each corresponding time points after transition to heterotrophic condition (autotrophic condition: SRR7650945 and SRR7650944; 3 hour of growth in heterotrophic culture: SRR7650947, SRR7650946; 6 hour of growth in heterotrophic culture: SRR7650952, SRR7650943; 12 hour of growth in heterotrophic culture: SRR765095, SRR7650950; 24 hour of growth in heterotrophic culture: SRR7650942, SRR7650941; 72 hour of growth in heterotrophic culture: SRR7650949, SRR7650948).

Processing of RNA-seq data

The quality control of raw data was performed using FastQC software version v.0.11.5. Low-quality reads were filtered out using Trimmomatic v.0.36 [15] with the following parameters [LEADING: 30, TRAILING: 3, SLIDINGWINDOW: 4:20, and MINLEN: 30]. Low-quality bases (quality scores <3) were trimmed out. Trimmed reads were mapped to the A. protothecoides genome (assembly ASM73321v1) using TopHat2 software version 2.0.12 [16] with following parameters [–read-mismatches 5–read-gap-length 3–read-edit-dist 5–library-type = fr-firststrand–splice-mismatches 0]. Then, sorted Alignment/Map (BAM) files were subjected to Picard software version 2.5.0 (https://broadinstitute.github.io/picard/) to remove potential duplicates.

Analysis of differentially expressed genes

HTSeq-count [17] was used to quantify gene expression at each sample. Prior study have corroborated that negative binomial based methods are more efficient to identify DEGs in experiments with less biological replicates [18]. Therefore, created count matrix was subjected to negative binomial based methods implemented in the Bioconductor DESeq2 package version 1.10.1 [19]. To avoid biases in fold change, low read counts (<5) across all samples were filtered out. Normalization of count data for all samples was performed using the default settings in DESeq2. Wald’s test was applied to determine the log2-fold change between AG and HG condition in five time points viz, 3h, 6h, 12h, 24h, and 72h. Significant differential expression was determined as a fold change ≥ |2| with a false discovery rate (FDR) corrected p-value cut off of ≤ 0.05 [20]. Heat map analysis was also performed to visualize of common DEGs at different time points grown in heterotrophic condition.

Genome-wide analysis of alternative splicing

Processed RNA-seq data were also used to genome wide study of alternative splicing changes of microalgae at heterotrophic condition. Due to reduced illative power, analysis of alternative splicing using the transcript assembly is challenging, therefore, differential exon usage as evidence for alternative splicing was analyzed by using the state-of-the-art software DEXSeq version 1.16.10 [17]. Exon expression was quantified by using a modified HTSeq-count phyton script that was provided in the DEXSeq package [17].

To test for differential expression usage in each gene across all the samples during the transition from autotrophic to heterotrophic condition following interaction model (~ sample + exon + condition:exon) and a likelihood ratio test were applied. Significant DEU was considered as a false discovery rate (FDR) corrected p-value ≤ 0.05 [20].

ASTALAVISTA program [21] was also applied to detect AS landscape, including Intron retention (IntronR), alternative 3' splice acceptor (AltA), alternative 5′ splice donor (AltD), Exon skipping (ExonS) and complex event. The output from landscape analysis can be found in the S1 and S2 Tables for heterotrophic and autotrophic condition, respectively.

Functional analysis

To examine the functional impacts of the DEG and DAS genes, GO analysis was performed to screen for the enriched GO terms with Fisher’s exact test using the Web Gene Ontology Annotation Plot (WEGO) (http://wego.genomics.org.cn/) [22]. Enriched GO terms were categorized in three Biological Process (BP), Molecular Function (MF), and Cellular Component (CC) group under p-value < 0.05.

Carotenoids biosynthesis‑related AS genes

To identify carotenoids biosynthesis‑related genes of A. protothecoides, corresponding genes from Chlamydomonas reinhardtii and Chlorella vulgaris were downloaded from the uniprot database (https://www.uniprot.org/). Then, BLASTP was used to align protein sequences of the A. protothecoides (available at https://www.ncbi.nlm.nih.gov/genome/) as reference and the downloaded sequences as the query with cutoff E value 1E–5.

Validation of alternative splicing

Additional independent A. protothecoides RNA-Seq dataset PRJNA289168 based on the same conditions as the training set was downloaded from ENA database and used as a test set to verify existence and patterns of splicing in candidate DAS genes at heterotrophic condition. This dataset included three autotrophic samples grown in Bristol's salts plus 0.1% (w/v) proteose peptone and nine samples grown in papaya-based heterotrophic cultures (Bristol's salts plus 12.5% (v/v) Rainbow papaya juice.

Results

Dynamic transcriptome changes at heterotrophic growth mode

In total, 3614, 2739, 2869, 3542, and 2234 DEGs were found on hour 3, hour 6, hour 12, hour 24, and hour 72, respectively (Fig 1A). Among the mentioned DEGs, 705 of which were shared in all five time points (Fig 1B). About 54% (hour 3: 45%, hour 6: 58%, hour 12: 47%, hour 24: 65%, hour 72.56%) DEGs were up regulated in HG cells (Fig 1C).

[Figure omitted. See PDF.]

Fig 1.

Profiles of DEGs in microalgae after transition to heterotrophic condition, frequency of up regulated and down regulated gene (A), overlap of DEGs (B) and expression pattern (C) of common DEGs at five time point after transition to heterotrophic condition.

https://doi.org/10.1371/journal.pone.0234710.g001

Top up- and down-regulated genes and their function annotations on hour 3, hour 6, hour 12, hour 24, and hour 72 of growth in heterotrophic condition are listed in S3 Table. Results showed that ammonium transporter 1 member 3 (AMT1-3) and purple acid phosphatase 15 (PAP15) are among the top up-regulated and superoxide dismutase 2 (SOD2) is among the top down-regulated genes.

Genome-wide analysis of alternative splicing in heterotrophic condition

In current study, besides differential expression analysis, genome-wide survey of AS process in A. protothecoides was also performed. 5229 AS events in autotrophic condition and 11618 AS events in heterotrophic condition were detected, which indicated that heterotrophic condition dramatically induce the AS process in A. protothecoides (S4 Table). The overall statistics of the shared/unique DAS at five time points in comparison with DEGs are shown in Fig 2. 660 DAS were common at different time points (Fig 2). Five types of AS events viz, IntronR, AltA, AltD, ExonS, and others, were considered in current study (Fig 2). Results showed that IntronR and AltD are the most prevalent type of AS in heterotrophic and autotrophic growth condition (Fig 2). On average, 28% of the AS events were IntronR (28%), followed by others (27%), AltD (24%), and AltA (15%), with ExonS (5%) being the least frequent type of AS event at heterotrophic growth mode (Fig 2). Validation of occurrence and pattern of AS process for candidate DAS genes were performed by using another RNA-Seq dataset in heterotrophic and autotrophic growth mode for A. protothecoides. AS events analysis in additional datasets detected the similar rate and pattern for AS underwent genes during transition from autotrophic to heterotrophic condition.

[Figure omitted. See PDF.]

Fig 2.

Profiles of DAS in microalgae after transition to heterotrophic condition, overlap of DAS genes at five time points (A), catalogue of AS event in autotrophic and heterotrophic condition (B) and overlap of DAS and DEGs (C) in heterotrophic condition.

https://doi.org/10.1371/journal.pone.0234710.g002

Functional impact of DEG and DAS genes

GO annotation showed that DEG and DAS gene sets had similar GO term coverage indicated that these two gene sets participated in the same biological process and cooperated with each other, producing the response to heterotrophic condition. Comparison of GO for DAS and DEG sets shows that the frequency of GO terms is different (Fig 3).

[Figure omitted. See PDF.]

Fig 3. Functional enrichment analysis of DEG and DAS genes.

Significantly (P-value<0.05) enriched gene ontology in three categories including cellular component, molecular condition and biological process are shown in below of the figure.

https://doi.org/10.1371/journal.pone.0234710.g003

Genes categorized by ‘binding’ and ‘response to stimuli’, for example, were higher in DAS than DEG set. Meanwhile, ‘oxidoreductase activity’, ‘transcription regulatory activity’, and ‘cofactor binding’ were highly enriched in DEGs. In molecular function category, binding (GO: 1901363), protein binding (GO: 0005515), catalytic activity (GO: 0003824), molecular function regulator (GO: 0098772), antioxidant activity (GO: 0016209), molecular transducer activity (GO: 0060089), molecular transducer activity (GO: 0060089) and transcription regulator activity (GO: 0140110) were significantly enriched in DAS genes (P-value < 0.05). In the GO term binding (GO: 0005488), heterocyclic compound binding (GO: 1901363) including nucleic acid binding (GO: 0003676), vitamin B6 binding (GO: 0070279) and tetrapyrrole binding (GO: 0046906), protein binding (GO: 0005515) including protein dimerization activity (GO: 0046983), unfolded protein binding (GO: 0051082), heat shock protein binding (GO: 0031072), transcription factor binding (GO: 0008134) were significant.

Alternative splicing and expression patterns of the carotenoids biosynthesis related genes

We analyzed the genes involved in carotenoid biosynthesis. Ten A. protothecoides genes that encode enzymes in the carotenoid biosynthetic pathways, including geranylgeranyl pyrophosphate synthase 1(GGPS1), phytoene synthase (PSY), phytoene desaturase (PDS), ζ-carotene desaturase (ZDS), lycopene epsilon cyclase (LCYE), lycopene beta cyclase (LCYB1), β-carotene hydroxylase (CHYB), P450 hydroxylases (P450-CHY), Carotenoid isomerase (CrtISO), and zeaxanthin epoxidase (ZEP) were selected for alternative splicing and expression analysis. The heat map of the carotenoid biosynthesis related genes was drawn according to fold change values (Fig 4).

[Figure omitted. See PDF.]

Fig 4. Expression profile of carotenoid biosynthesis related genes after transition to heterotrophic condition.

https://doi.org/10.1371/journal.pone.0234710.g004

Except the P450b-CHY, ZDS, and LCYB1, expression of the most carotenoid biosynthetic pathway genes decreased at heterotrophic growth. We next tried to investigate whether heterotrophic growth condition would induce alternative splicing of carotenoids biosynthesis related genes. Results indicated that heterotrophic condition affects AS pattern of carotenoids biosynthesis related genes. Comparison of AS pattern in autotrophic and heterotrophic growth condition showed that PDS, LCYB1, CrtISO, ZEP, and LCYE undergo AS process at autotrophic condition, whereas in heterotrophic condition, AS events were detected for the PDS, P450-CHY, ZDS, and CrtISO, highlighting the regulatory impacts of growth condition changes on splicing pattern of carotenoids biosynthesis pathway.

Discussions

Recently developed high throughput sequencing technology so called Next Generation Sequencing (NGS) has allowed effective genome-wide detection of alternative splicing events. In current study, we used RNA-seq data to uncover features of alternative splicing and gene expression in A.protothecoides grown in heterotrophic condition. Comparison of DEG and DAS showed that it was relatively small overlap between DAS and DEG genes (84 genes) indicated that AS accomplish as a distinct adaptive regulatory layer in heterotrophic condition. Functional enrichment analysis revealed that the DAS genes were involved mainly in several biological processes, including cellular metabolic process, nitrogen compound metabolic process, and biosynthetic process. Further analysis revealed that DAS genes were involved tetrapyrrole binding (GO: 0046906). Among different genes involved in tetrapyrrole binding, the gene encodes the multidrug resistance (MDR)-like ABC transporters undergoes AS events in heterotrophic condition. This protein has been shown to catalyze the primary active export of auxin [23] highlighting the functional importance of AS event in signaling rewiring during growth mode change. Among different stress related genes, DnaJ-like proteins belong to heat shock protein 40 (Hsp40) families were undergoing AS events in heterotrophic condition. It has been demonstrated that DnaJ-like proteins participate in regulation of protein folding, protein transport and cellular responses to stress [24]. Intriguingly, transcription regulator activity significantly enriched in DAS genes. More dissection of DAS genes that were categorized in transcription regulator activity, showed that DNA-binding transcription factor activity, transcription initiation factor activity are mostly enriched in DAS genes, whereas these categories were not found in DEGs which indicated that AS can serve as an independent mechanism for gene regulation in response to growth mode change [25, 26].

The most striking result to emerge from the data is that growth in heterotrophic condition affects AS pattern of carotenoids biosynthesis pathway. Among different genes contributed in carotenoids biosynthesis pathway, PDS, P450-CHY, ZDS and CrtISO underwent AS event. AltD, AltA, IntronR, and AltD events were detected for PDS, P450-CHY, ZDS, CrtISO, respectively. More interestingly, similar pattern of AS were also observed in training data set for carotenoids biosynthesis related genes.

PDS, a key enzyme in carotenoid biosynthesis, are needed in the conversion from phytoene to lycopene [27]. It has 8 exons (E) of which 3 exons show significant DEU and 5 exons are conserved usage at heterotrophic condition in compared with autotrophic condition. E3 and E8 of PDS exhibit time point-specific exon usage at heterotrophic condition. Only E1 and E2 in this gene were consistently showed DEU at all-time points (Fig 5). In agreement with our finding, it has been previously suggested that alternative transcripts could regulate the cellular concentrations of PDS, depending on the physiological or environmental conditions [28].

[Figure omitted. See PDF.]

Fig 5. Gene models illustrating AS patterns of PDS gene.

Exon usage was used to measure exon expression after normalizing for total gene expression biases. Significant differentially expressed exons are colored in magenta (p.adjust<0.05) in 6h (A), 12h (B), 24h (C) and 48h (D) after transition to heterotrophic condition.

https://doi.org/10.1371/journal.pone.0234710.g005

ZDS and CrtISO convert ζcarotene to lycopene, which is the substrate for two competing lycopene cyclases, ε-LCY and β-LCY [29]. It has been reported that overexpression of ZDS significantly increase the lutein content and salt tolerance [30]. Finding of our study showed that ZDS undergoes AS event (IntronR) at heterotrophic growth condition (time point 6, 12, and 72 h). ZDS has ten exons of which three exons showed significant DEU at heterotrophic condition. Exon 4 and 7 consistently showed DEU at hour 6, 12, and 72, after transition to the heterotrophic condition, however, DEU of exon 10 was only found at hour 72.

Regarding the CrtISO, the number of exons with DEU was varied from hour 3 to 72. It is apparent from the results that at 3 and 6 hour after transition to the heterotrophic condition, E 5 and E 11 show DEU, it is whilst at hour 12 and 72 the number of exons with DEU increased to 4 (E1, E2, E7, and E 11) and 5 (E1, E2, E5, E7 and E11), respectively. Whereas no significant expression was identified at hours 24 and 72, highlighting the importance of post-transcriptional modification processes on carotenoids biosynthesis pathway regulation during adaptation to different growth condition. In line with our hypothesis, it has been suggested that the translational control mechanisms is interpreted in terms of flux adjustments needed in response to retrograde signals stemming from intermediates of the plastid-localized carotenoid biosynthesis pathway [31].

Conclusions

In current study, we present a RNA-seq analysis of a comparison of the transcriptional and post-transcriptional responses of A. protothecoides grown under autotrophic and heterotrophic conditions. More AS events in the heterotrophic condition were identified than in the autotrophic condition. There is relatively small overlap between DAS and DEG gene indicated that AS accomplish as a distinct adaptive regulatory layer in responses to growth condition change. To the best of our knowledge, this is the first study in which alternative splicing has investigated in the microalgae grown on heterotrophic and autotrophic conditions. The results of current study provides new insight into underlying regulatory responses of microalgae grown in modified condition to optimized production of high value metabolites at industrial scale.

Supporting information

[Figure omitted. See PDF.]

S1 Table. Alternative splicing landscape analysis in autotrophic condition.

https://doi.org/10.1371/journal.pone.0234710.s001

(TXT)

S2 Table. Alternative splicing landscape analysis in heterotrophic condition.

https://doi.org/10.1371/journal.pone.0234710.s002

(GTF)

S3 Table. Differentially expressed genes and their corresponding fold changes at 6, 12, 24, and 72 hour after transition to heterotrophic condition growth in heterotrophic condition.

https://doi.org/10.1371/journal.pone.0234710.s003

(XLSX)

S4 Table. Differentially spliced genes at 6, 12, 24, and 72 hour after transition to heterotrophic condition growth in heterotrophic condition.

https://doi.org/10.1371/journal.pone.0234710.s004

(XLSX)

Acknowledgments

We would like to thank Dr Asgar Panahi for reading of the manuscript.

Citation: Panahi B, Hejazi MA (2020) Integrative analysis of gene expression and alternative splicing in microalgae grown under heterotrophic condition. PLoS ONE 15(6): e0234710. https://doi.org/10.1371/journal.pone.0234710

References

1. Harun R, Singh M, Forde GM, Danquah MK (2010). Bioprocess engineering of microalgae to produce a variety of consumer products. Renewable and Sustainable Energy Reviews. 14(3):1037–47.

2. Gladue RM, Maxey JE (1994). Microalgal feeds for aquaculture. Journal of Applied Phycology. 6(2):131–41.

3. Akubude VC, Nwaigwe KN, Dintwa E (2019). Production of biodiesel from microalgae via nanocatalyzed transesterification process: A review. Materials Science for Energy Technologies. 2(2):216–25.

4. Xiao Y, He X, Ma Q, Lu Y, Bai F, Dai J, et al (2018). Photosynthetic accumulation of lutein in Auxenochlorella protothecoides after heterotrophic growth. Marine drugs. 16(8):283–292.

5. Klamczynska B, Mooney WD (2017). Heterotrophic microalgae: a scalable and sustainable protein source. In Sustainable protein sources. Jan 1 (pp. 327–339). Academic Press.

6. Panahi B, Farhadian M, Hejazi MA (2020). Systems biology approach identifies functional modules and regulatory hubs related to secondary metabolites accumulation after transition from autotrophic to heterotrophic growth condition in microalgae. Plos one. 15(2):e0225677. pmid:32084664

7. Panahi B, Mohammadi SA, Khaksefidi RE, Mehrabadi JF, Ebrahimie E (2015). Genome-wide analysis of alternative splicing events in Hordeum vulgare: Highlighting retention of intron-based splicing and its possible function through network analysis. FEBS letters. 589(23):3564–3575. pmid:26454178

8. Filichkin SA, Priest HD, Givan SA, Shen R, Bryant DW, Fox SE, et al (2010). Genome-wide mapping of alternative splicing in Arabidopsis thaliana. Genome research. 20(1):45–58. pmid:19858364

9. Beisel NS, Noble J, Barbazuk WB, Paul AL, Ferl RJ 2019. Spaceflight-induced alternative splicing during seedling development in Arabidopsis thaliana. NPJ microgravity. 5(1):1–5.

10. Ling Y, Alshareef S, Butt H, Lozano‐Juste J, Li L, Galal AA, et al (2017). Pre‐mRNA splicing repression triggers abiotic stress signaling in plants. The Plant Journal. 89(2):291–309. pmid:27664942

11. James AB, Syed NH, Bordage S, Marshall J, Nimmo GA, Jenkins GI, et al (2012). Alternative splicing mediates responses of the Arabidopsis circadian clock to temperature changes. The Plant Cell. 24(3):961–81. pmid:22408072

12. Kianianmomeni A, Ong CS, Rätsch G, Hallmann A (2014) Genome-wide analysis of alternative splicing in Volvox carteri. BMC genomics. 15(1):1117–1122.

13. Shang X, Cao Y, Ma L (2017). Alternative splicing in plant genes: a means of regulating the environmental fitness of plants. International journal of molecular sciences. 18(2):432–442.

14. Huang B, Xu Y (2015). Cellular and molecular mechanisms for elevated CO2–regulation of plant growth and stress adaptation. Crop Science. 55(4):1405–1424.

15. Bolger AM, Lohse M, Usadel B (2014). Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. 30(15):2114–20. pmid:24695404

16. Kim D, Pertea G, Trapnell C, Pimentel H, Kelley R, Salzberg SL (2015). TopHat2: accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions. Genome biology. 14(4): 36–40.

17. Anders S, Pyl PT, Huber W (2015). HTSeq—a Python framework to work with high-throughput sequencing data. Bioinformatics. 31(2):166–169. pmid:25260700

18. Schurch NJ, Schofield P, Gierliński M, Cole C, Sherstnev A, Singh V, et al (2016). How many biological replicates are needed in an RNA-seq experiment and which differential expression tool should you use?. Rna. 22(6):839–51. pmid:27022035

19. Love MI, Huber W, Anders S (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome biology. 15(12):550–562. pmid:25516281

20. Benjamini Y, Hochberg Y (1995). Controlling the false discovery rate: a practical and powerful approach to multiple testing. Journal of the Royal statistical society: series B (Methodological). 57(1):289–300.

21. Foissac S, Sammeth M (2007). ASTALAVISTA: dynamic and flexible analysis of alternative splicing events in custom gene datasets. Nucleic acids research. 35:297–299.

22. Ye J, Fang L, Zheng H, Zhang Y, Chen J, Zhang Z, et al (2006). WEGO: a web tool for plotting GO annotations. Nucleic acids research. 34(suppl_2):W293–297.

23. Santelia D, Vincenzetti V, Azzarello E, Bovet L, Fukao Y, Düchtig P, et al (2005). MDR‐like ABC transporter AtPGP4 is involved in auxin‐mediated lateral root and root hair development. FEBS letters. 579(24):5399–53406. pmid:16198350

24. Liu J, Li X, Tang X, Zhou B (2016). Identification and characterization of a DnaJ gene from red alga Pyropia yezoensis (Bangiales, Rhodophyta). Chinese journal of oceanology and limnology. 34(2):405–411.

25. Panahi B, Abbaszadeh B, Taghizadeghan M, Ebrahimie E (2014). Genome-wide survey of alternative splicing in Sorghum bicolor. Physiology and Molecular Biology of Plants. 20(3):323–9. pmid:25049459

26. Pimentel H, Parra M, Gee SL, Mohandas N, Pachter L, Conboy JG (2017). A dynamic intron retention program enriched in RNA processing genes regulates gene expression during terminal erythropoiesis. Nucleic acids research. 44(2):838–851.

27. Huang J, Liu J, Li Y, Chen F (2008). Isolation and characterization of the phytoene desaturase gene as a potential selective marker for genetic engineering of the astaxanthin-producing green alga Chlorella zofingiensis (chlorophyta)(1). Journal of Phycology. 44(3): 684–690. pmid:27041426

28. Lodato P, Alcaino J, Barahona S, Retamales P, Cifuentes V (2003). Alternative splicing of transcripts from crtI and crtYB genes of Xanthophyllomyces dendrorhous. Applied Environmental Microbiology. 69(8):4676–4682. pmid:12902257

29. Cazzonelli CI, Pogson BJ (2010). Source to sink: regulation of carotenoid biosynthesis in plants. Trends in plant science. 15(5):266–274. pmid:20303820

30. Li R, Kang C, Song X, Yu L, Liu D, He S, et al (2017). A ζ-carotene desaturase gene, IbZDS, increases β-carotene and lutein contents and enhances salt tolerance in transgenic sweetpotato. Plant Science. 262:39–51. pmid:28716419

31. Álvarez D, Voß B, Maass D, Wüst F, Schaub P, Beyer P, et al (2016). Carotenogenesis is regulated by 5′ UTR-mediated translation of phytoene synthase splice variants. Plant physiology. 172(4):2314–2326. pmid:27729470