Transcriptomes are highly sensitive, with widespread variability among experimental units exposed to the same treatment detected in even highly controlled environments (Cortijo et al., 2019). Minimizing undesired transcriptional variation can maximize the success of RNA-seq experiments.

Preprocessing comprises the first step for any sequencing analysis. This stage involves the removal of technical artifacts such as adaptors, PhiX sequence, rRNA sequences, assessing sequence quality and if necessary, quality trimming. Commonly, sequencing centers perform the preprocessing steps of the RNA-seq pipeline, but it is best practice to perform a data quality check in-house. A multitude of tools are available for the preprocessing steps. Quality of raw sequences can be checked via tools such as FastQC (Andrews, 2010) and MultiQC (Ewels et al., 2016) or fastp (Chen S. et al., 2018). Preprocessing tools such as such as Cutadapt (Martin, 2011), Trimmomatic (Bolger et al., 2014), BBTools (Bushnell, 2022), and fastp (Chen et al., 2018) can be used for contamination removal or quality trimming. rRNA contamination can be identified and removed with tools like BBTools and SortMeRNA (Kopylova et al., 2012).

Preprocessing steps should be implemented with caution, as they impact downstream analyses. A more stringent read trimming shortens reads and thus influences mapping and the estimates of expression levels for genes and their isoforms, ultimately impacting the differential expression analysis (Williams et al., 2016). For de novo transcriptome analysis, a stringent trimming strategy has also been observed to produce incomplete transcriptome reconstruction (MacManes, 2014; Mbandi et al., 2014). To ensure the greatest amount of information is retained for analysis, employing less stringent read-trimming, or no read-trimming at all is suggested (Del Fabbro et al., 2013; MacManes, 2014; Mbandi et al., 2014; Williams et al., 2016).

The next step in the RNA-seq analysis pipeline is to align reads using a reference genome or de novo methods. Similar to preprocessing steps, there are many bioinformatic tools available for aligning reads, and they are categorized by whether they require a reference genome or a de novo transcriptome. In this section we provide a brief overview of the different methods available for both reference and de novo-based approaches, with a focus on the types of biological questions each method can address and the limitations of each approach.

Alignment-based reconstruction occurs when a reference genome is available. The approach is very similar to aligning genomic reads, but splicing events need to be considered for RNA-seq data (Mehmood et al., 2020). Some of the software available to perform the reference-based alignment includes HISAT2 (Kim et al., 2019) and STAR (Dobin et al., 2013). Reference-based alignment is split into two parts: reference genome indexing and alignment of reads to the indexed reference genome. Additionally, some alignment tools can incorporate the discovery of novel exon-exon junctions like HISAT2 (Kim et al., 2019) and RsubRead (Liao et al., 2019).

One drawback to reference-based alignments is reference-bias: if there are sequences present in the RNA-seq data that are not present in the reference data, the data will not align and will be lost for downstream analyses unless alternative mining is performed (Yang et al., 2022). A reference-based transcript reconstruction can be performed post-alignment by software like Cufflinks (Trapnell et al., 2010), StringTie2 (Kovaka et al., 2019) and Bookend (Schon et al., 2022). Alternatively, new approaches that align reads to pangenomes are available. Pangenomes store population-level genetic variation into a graph-based structure rather than a single linear genome (Eizenga et al., 2020), allowing for improved read alignment, including for haplotype-aware RNA-seq read alignment (Sibbesen et al., 2023).

When a reference genome or pangenome is not available for your species of interest, researchers can choose to perform a de novo assembly transcriptome. De novo transcriptome assembly can be performed by tools as Trinity (Grabherr et al., 2011) and TransLiG (Liu et al., 2019). The same tools for a reference genome-based approach are used to align RNA-Seq to a de novo transcriptome with minor modifications – e.g., not using the splice-aware function of the aligners.

The next step in the bioinformatics pipeline is counting reads mapped. For reference-based approaches tools such as HT-Seq (Anders et al., 2015) or featureCounts can be used (Liao et al., 2014). Also, transcript quantification using a reference transcriptome can be performed using alignment-free methods like Kallisto (Bray et al., 2016) and Salmon (Patro et al., 2017). It is worthwhile mentioning that RNA-Seq strandedness impacts quantification, identification of isoforms and de novo transcriptome assemblies. The concept refers to the strategy employed in library construction, wherein stranded library preparations maintain the transcript directionality. Hence, researchers should be aware of the kit utilized for library construction when performed the aforementioned steps of RNA-seq analysis. When the information of strandedness and direction of strandedness is not known, tools like how_are_we_stranded_here (Signal and Kahlke, 2022) can help determine the strandedness of paired-end libraries. Finally, tools like tximport (Soneson et al., 2015) and tximeta (Love et al., 2020) can be used to summarize the quantification of transcript abundances into an expression matrix. For a detailed review about alignment and quantification, we refer the reader to Van den Berge et al. (2019).

Measuring plant phenotypic plasticity in extreme environments can aid in knowledge to develop plants for future abiotic and biotic environmental conditions arriving with climate change (Gage et al., 2017; Monforte, 2020). Using various molecular mechanisms, plants respond to changing environmental conditions by altering their physiology and development (Lachowiec et al., 2015). Molecular plant plasticity enables adaptation to climatic shifts and predicts an individual’s survival success (Nicotra and Davidson, 2010; Nicotra et al., 2010; Fox et al., 2019; Anderson and Song, 2020; Pazzaglia et al., 2021), warranting further investigation (Brooker et al., 2022).

Increasingly studies use RNA-seq to understand molecular plant plasticity (Kumar et al., 2022; Sreeratree et al., 2022). The primary approach analyzes differential expression (DE), which evaluates the transcriptional abundance across conditions through simultaneous statistical testing for significant changes in expression levels in all detected genes, transcripts, or different usage of transcripts/exons (Soneson et al., 2015; Van den Berge et al., 2019). Software like edgeR (Robinson et al., 2010), DESeq2 (Love et al., 2014), and limma (Ritchie et al., 2015; Stark et al., 2019) all provide robust DE analyses. For reviews of the main aspects in the differential expression analysis we refer the reader to Costa-Silva et al. (2021); Stark et al. (2019) and Van den Berge et al. (2019).

Differential expression analysis incorporates a matrix of features’ expression levels and knowledge about the experimental design. Tests for identifying differentially expressed genes (DEGs) rely on contrasting conditions, such as different tissues, genotypes, and conditions. Exploring the up- or downregulation of genes under a stress condition relative to a control condition indicates how a plant combats a stressor and how the stressor harms the plant. For example, in the sorghum lateral root apex, low levels of phosphorus caused major expression changes in the lateral root apices, which correlated with enhanced lateral root growth. Specifically, the low-phosphorus-induced genes encoded proteins with functions in nutrient responses and contribution to phosphorus metabolism (Gladman et al., 2022). Contrasting genotypes with different performances under stress enables identifying potential mechanisms of stress resistance (Yue et al., 2016; Zhao et al., 2021).

Allele-specific expression (ASE) describes the phenomenon where alleles within a particular genetic feature (e.g. gene, transcript) have significant differences in their expression levels (Castel et al., 2015). The expression of alleles can be assessed by comparing the expression of genes of a certain genotype with its parents or by identifying polymorphisms and quantifying the expression of each allele. ASE is especially informative in understanding hybrid crops where different parental genotypes are combined (Bell et al., 2013; Shao et al., 2019). ASE analysis also reveals the processes of genetic imprinting, tissue- and stress-specific alleles, as well as the evolution of species.

In wheat seeds, expression of homeologs and alleles is differentially controlled and consideration of each copy of a gene is relevant. In the endosperm, genes exhibited subgenome dominance in particular functions (Pfeifer et al., 2014). Further, imprinted genes were identified more frequently in developing endosperm relative to other tissues, and imprinted gene expression patterns were conserved through wheat evolution (Yang et al., 2018). This imbalanced expression of maternal and paternal alleles and subgenome dominance supports proper seed development.

ASE analysis uncovered genes that were targets of selection during domestication with implications for plant sciences. Lemmon et al. (2014) identified that maize and teosinte diverged in gene expression especially due to cis regulation, where the expression of maize alleles is favored in F1 hybrids between the species. Genes with cis and cis plus trans divergent regulation were enriched among putative targets of selection. In rice, a F1 hybrid of genotypes representing two major subpopulations exhibit enrichment of genes ASE in genomic regions of signatures for domestication or artificial selection (Shao et al., 2019). Moreover, the limited ASE in sugarcane internodes predominates as genotype-specific phenomenon, favoring high dosage alleles and purging the expression of potentially deleterious alleles (Margarido et al., 2022). Detecting ASE and regulatory mechanisms of ASE can inform our understanding of plant evolutionary history.

Alternative splicing of precursor mRNAs leads to diversification of the functions of a single gene. Transcriptional or isoform switching refers to a shift in the presence or dominance of transcripts in different samples, including across cell types, development, genotypes, and environments. Within plants, the most common form of alternative splicing is intron retention (Chamala et al., 2015), in contrast to animals where exon skipping is most detected (Kim et al., 2007). When the chosen alignment and counting approaches (see Section 3.2) enable distinguishing transcripts, examination of transcript switching is possible by utilizing the differential expression of individual transcripts or using accessible tools for assessing transcriptional switching (Qiu et al., 2021).

Genome-wide surveys of alternative splicing demonstrate the potential impact of alternative splicing events, with tens of thousands of alternative splice forms detected that are evolutionarily dynamic across angiosperms (Chamala et al., 2015). Across a population of over 350 inbred maize lines, variation in alternative splicing was detected, highlighting that connecting genotype to phenotype can be better informed by considering the expression of particular splice forms (Chen Q. et al., 2018).

Co-expression summarizes large-scale transcriptomics to infergene regulatory networks by identifying modules of genes with similar expression patterns across multiple samples. Co-expression analyses can suggest target genes of interest and corroborate GWAS results and expand potential genetic markers from those findings as well. Putative functions then can be assigned to non-annotated genes if the majority of the genes in a module share similar biological functions, the guilt-by-association principle (Serin et al., 2016; Rao and Dixon, 2019). Co-expression network analyses split into two main approaches: (i) non-targeted—a network based on the topological structure using the relationship of all pairs of genes or (ii) targeted—the use of bait genes as prior information to define network connections. Varied inputs are used to build networks, including replicates of multiple treatments, averaging the expression of replicates grown for a treatment, or defining networks separately to single levels of the experimental factor (Cortijo et al., 2020).

A fundamental step for inferring the co-expression of genes from large-scale transcriptomic data is the use of similarity measures, including correlation and mutual information methods (Ma and Wang, 2012; Huang et al., 2017). Correlation methods selection is based on data type, and common coefficients include Pearson’s, Spearman’s and Gini’s (Huang et al., 2017). Gini’s correlation, for example, is advantageous for nonnormally distributed RNA-Seq data, robust against outliers and small sample sizes (Ma and Wang, 2012). A common pipeline to construct a co-expression network involves the calculation of a similarity matrix that is then filtered based on a threshold to select gene pairs; then an adjacency matrix can be calculated and subsequently a clustering algorithm is used to group genes into modules (Serin et al., 2016).

Studies use co-expression networks to find hubs—genes with high network connectivity—and understand their role in biological pathways. For example, sugarcane hub genes changed across four stages of development in the networks of 10-month-old compared to 6-month-old apical culms (Hosaka et al., 2021), revealing candidates related to cell wall and stress and three transcription factors (TFs) potentially acting as regulators of those processes. Co-expression networks may also uncover functions for uncharacterized proteins. De Vega and colleagues (De Vega et al., 2021) inferred TF targets in Miscanthus hybrids that were enriched with carbohydrate metabolism, secondary metabolism, and the generation of precursor metabolites. They also found two TF that linked a core and a loop subnetwork, the last composed mostly by TFs linked to uncharacterized genes. Thus, regulatory co-expression networks are useful tools to identify targets of TFs, which can be important targets for biotechnology and propose functions for poorly understood proteins (Simons et al., 2006; Fröschel et al., 2019).

Time-series expression data impose a challenge for understanding the dynamics of the coregulation of genes. The complex relationships arising due to time can be detected by gene co-expression measures that account for local dependence structures in the expression patterns (Wang et al., 2014). The dynamic network biomarker approach (Chen et al., 2012) aims to find a subnetwork of strongly correlated genes just before a critical transition – identified as a tipping point. With this approach researchers found genes at the tipping point for response to stress or ripening of fruits (Wang and Zhang, 2021; Wang et al., 2022).

Co-expression networks can also provide additional data for determining genes involved in regulatory networks. Co-expression combined with ChIP-seq can also identify the targets of a TF. Cortijo and colleagues (Cortijo et al., 2020) identified novel regulatory targets in the Arabidopsis thaliana circadian clock by combining modules with genes co-expressed with known TFs and ChIP-seq data. Targets of PSEUDO-RESPONSE REGULATOR 5—a core component of the circadian clock—were found in a module showing enrichment for photosynthesis. For a complete review of co-expression networks in plant biology, the reader is referred to (Rao and Dixon, 2019) and (Serin et al., 2016).

Enrichment analysis tests if any functional group is over- or underrepresented by a list of genes of interest – e.g., DEGs, co-expressed genes or ASE genes. Enrichment analysis of genes identified by differential expression or modules in co-expression networks are useful to understand if genes in the module are related to similar functions. It can suggest that non-annotated genes likely participate in the same biological pathways as known genes, which can lead to the identity of causal genes (Serin et al., 2016). Functional pathways are represented by ontologies in different frameworks like Gene Ontology (GO) categories (Ashburner et al., 2000), Kyoto Encyclopedia of Genes and Genomes (KEGG) orthology (KO) (Kanehisa et al., 2017), and MapMan4 bin categories (Schwacke et al., 2019), which all capture varied functions of genes.

For plants with reference genomes, functional annotation is provided by databases along with the nucleotide sequence and the structural annotation. For de novo assemblies this procedure requires comparisons with nucleotide or protein databases, retrieving the functional categories from the hits and associating them to genes/transcripts of the reference. While functional annotation of a reference – genome or transcriptome – in terms of GOs, Enzyme Codes (ECs) or KOs can be performed by tools like OmicsBox (BioBam Bioinformatics, 2019) or Trinotate (Bryant et al., 2017), Mercator4 provides the specific annotation for MapMan (Schwacke et al., 2019). The KO/EC annotation can be linked to molecular networks by KEGG Mapper (Kanehisa and Sato, 2020) to visual representation of the pathways. MapMan4 also has its own visualization of functional pathways where besides the category mapping file, fold change or expression values can be provided as input.

Finally, tests for enrichment of categories—over or underrepresentation—frequently make use of Fisher’s Exact test, Hypergeometrical test, permutation, or variation of those methods. There are many tools for functional enrichment analysis: goseq (Young et al., 2010), GSEA (Subramanian et al., 2005), OmicsBox (BioBam Bioinformatics, 2019), topGO (Alexa et al., 2006), clusterProfiler (Yu et al., 2012), and agriGO (Du et al., 2010).

Viruses are the most abundant biological entities on earth and are strongly associated with all organisms (Rosani and Gerdol, 2017; Mushegian, 2020). It is very common for viruses to contaminate tissue when preparing samples for sequencing. One meta-analysis found that over half of the 700 sequencing libraries examined had viral contamination (Asplund et al., 2019). RNA-seq is susceptible to viral contamination. However, most viral contamination is filtered out during bioinformatic processing.

Viral contamination in RNA-seq studies can provide informative results as well. Kamitani et al. (2016) looked at both plant-virus and virus-virus interactions in natural environments using RNA-seq. Other studies have used RNA-seq to screen for viral pathogens (Selitsky et al., 2020) and novel viral genome discovery (Rosani and Gerdol, 2017). However, some virus types require alternative rRNA depletion methods to be detected, such as non-polyadenylated genomes (Nagano et al., 2015). To capture the full range of viral genome varieties Nagano et al. (2015) optimized rRNA depletion methods for mRNA focused RNA-seq to detect varying types of viral signatures from natural environments. Viral contamination in transcriptomic data may provide insight to unknown biotic stressors or uncover novel defense mechanisms related to a previously unknown viral pathogen. Accessible software for identifying the species of origin allows for examining unmapped sequencing reads and examining expression patterns for additional species (Chen et al., 2020).

One feature of plant genomes that can influence the use of RNA-seq is the prevalence of polyploidy (Fu et al., 2016). A study of 203 modern crops identified that 17% have undergone polyploidization (Meyer et al., 2012). Of these, wheat and sugarcane are among the most cultivated crops globally (Weeks, 2017). The pairing of polyploid chromosomes during recombination varies along a spectrum from allopolyploidy to autopolyploidy, along with the presence of aneuploid chromosomes. Within allopolyploids, often formed through hybridization, chromosome copies behave in a diploid fashion and pair and segregate corresponding to species of origin, within subgenomes (Edger et al., 2018; Kuo et al., 2020). In contrast, homologous autopolyploid chromosomes pair at random, even between subgenomes (Spoelhof et al., 2017). The type of polyploidy can also vary along the length of a single chromosome. Polyploidization impacts the bioinformatic pipelines used to assess RNA-seq data, both at the level of processing reads and subsequent analyses. The presence of polyploidy can be problematic for RNA-seq analysis, from both a computational and biological standpoint.

Polyploidy increases computation complexity due to an increase in genome size, relative to a diploid genome and the increased repeat content. The increase in genome size, resulting from the presence of two or more subgenomes, means more memory is required to index and store the reference genome or transcriptome prior to read alignment. As such, analyzing large, polyploid genomes, such as wheat, may be prohibitive to those without expansive compute resources.

The degree of similarity among polyploid subgenomes is a driving force for how RNA-seq analyses may need to proceed differently compared to diploid genomes (Voshall and Moriyama, 2020). In the cases of subgenomes with low levels of divergence at the nucleotide level, little functional divergence may be expected across homeologs, the genes homologous to one another on each subgenome (Wang et al., 2017; Sigel et al., 2019). However, the phenomenon of subgenome dominance, or the tendency for genes to be expressed from a particular subgenomes, is commonly observed across plants (Schnable et al., 2011; Wang X. et al., 2016; Khan et al., 2020). For example, nonbalanced expression was identified in approximately ~30% of wheat homeologs, primarily through suppression of a single homeolog (Ramírez-González et al., 2018). Isolating the expression of homeologous genes may be particularly useful when using RNA-seq to complement QTL mapping studies (Yang et al., 2014). Therefore, distinguishing homeolog expression enables specific studies.

High levels of repeat content are also problematic from a computational viewpoint to interpret RNA-seq data. Repeat content can refer to both biological repeats (e.g., transposable elements, short sequence repeats) or genomic regions shared between subgenomes. Significant portions of widely grown cereals are comprised of repetitive sequences, excluding homeologs and duplicated genes: wheat-85% (Wicker et al., 2018), maize-85% (Schnable et al., 2009), barley-80% (Wicker et al., 2017), and rice-41% (Li et al., 2021). In the cases when these regions are expressed (Cavrak et al., 2014), they create issues computationally as there may be no way to confidently align reads to regions that are present more than once, resulting in a high number of multi-mapping reads. Multi-mapping reads are discarded by default by popular read counting programs. As such, the RNA-seq analysis of polyploid genomes may be hampered by the loss of a significant portion of data without explicit intervention.

Data loss due to repetitive sequences is problematic as repetitive elements can have important biological functions. For example, screening A. thaliana T-DNA insertion lines at the locations of transposable elements for seedling morphological responses to stresses uncovered functional roles of for over 90% of those tested (Joly-Lopez et al., 2017). Outside of transposable elements, short tandem repeats contribute to the repetitive content of the genome with functional consequences. The length of the short tandem repeat encoding a polyglutamine span in the protein EARLY FLOWERING 3 influenced flowering time across A. thaliana accessions (Undurraga et al., 2012). Inability to assign reads to these repeats means that functional polymorphisms are missing from analyses. Longer read technologies will improve reference genomes (Jiao et al., 2017; Michael et al., 2018; Kamal et al., 2022) and transcriptomes to enable study (Wang B. et al., 2016) of currently poorly characterized sequences.