Search
Submit Manuscript
1.6CiteScore
Advanced Search
MAP  >  Grass Research
2021 Volume 1
Article Contents
Abstract
Introduction
Materials and methods
Plant materials and gene library sequencing
De novo transcriptome assembly and quality analysis
Functional annotation of protein-coding transcripts
Identification of long non-coding transcripts
Identification of tentative reference genes
Evaluation of reliable reference genes
Quantification of cuticle biosynthesis-related transcripts
Statistical analysis
Results
Short reads and assembly quality
Homology searches
Functional categorization
Identification of lncrna
Identification of tentative reference genes
Expression stability of tentative reference genes
Identification of cuticle biosynthesis-related transcripts
Evaluation of reliable reference genes and quantification of cuticle biosynthesis-related transcripts
Discussion
S. thurberi transcriptome quality is similar to that reported for de novo assembly of non-model plants
S. thurberi transcriptome shows higher homology to sequences from species of the same family
Functional categorization analysis suggests a response to stress and an active cuticle biosynthesis in fruit pitaya peel
Lncrna from pitaya are similar in length and expression to that reported from other plants
Previously reported housekeeping genes were found to have the most stable expression during pitaya fruit development
Cutin biosynthesis could have a relevant role in the first stages of pitaya fruit development
Stabcg11 could be playing a relevant role in the transport of cuticle components during sweet pitaya fruit ripening
Conclusions
Author contributions
Data availability
Acknowledgments
Conflict of interest
Supplementary information
Rights and permissions
References
About this article
Next Previous
ARTICLE   Open Access    

Comparative transcriptomics and gene network analysis revealed secondary metabolism as preeminent metabolic pathways for heat tolerance in hard fescue

  • Department of Plant Biology, Rutgers, The State University of New Jersey, New Brunswick, NJ 08901, USA

More Information
  • The molecular mechanisms underlying genetic variations in heat tolerance, one of the important turfgrass traits, for fine fescue are not well-understood. In the present study, our objective was to identify molecular constituents and metabolic interactions involved in heat tolerance in two genotypes of hard fescue (Festuca trachyphylla) contrasting in heat tolerance by comparative transcriptomics and gene comparison network analysis. Two cultivars of hard fescue, 'Reliant IV' (heat-tolerant), and 'Predator' (heat-sensitive), were subjected to heat stress temperature at 35/30 °C (day/night) or maintained at optimal temperature at 22/18 °C (day/night) (non-stress control) for 21 d. At 14 and 21 d of heat stress, 'Reliant IV' maintained significantly higher photochemical efficiency (Fv/Fm), chlorophyll (Chl) content, and lower cell membrane electrolyte leakage (EL) compared to 'Predator', suggesting its superiority in heat tolerance. Comparative transcriptomic profiles, gene functional enrichment analysis, and weighted gene comparison network analysis revealed central hub genes (BBE22 and ALPLD) and their connecting genes involved in secondary metabolism for biosynthesis of oxylipins (LOX1 and LOX3), phenolic compounds (PAL2), and dhurrin (C79A1 and C71E1). These genes were up-regulated in heat-tolerant 'Reliant IV' under heat stress but not in heat-sensitive 'Predator', while a majority of heat-regulated genes involved in primary metabolism responded similarly to heat stress in both cultivars. Those unique genes in the secondary metabolic pathways enriched in only the heat-tolerant cultivar could be critical for mediating the protection of hard fescue against heat stress and are potentially useful as candidate genes or molecular markers for augmenting heat tolerance in other temperate species of grass.

  • Introduction

    Columnar cacti are plants of the Cactaceae family distributed across arid and semi-arid regions of America, with ecological, economic, and cultural value[1]. One trait that makes it possible for the columnar cactus to survive in the desert ecosystem is its thick epidermis covered by a hydrophobic cuticle, which limits water loss in dry conditions[1]. The cuticle is the external layer that covers the non-woody aerial organs of land plants. The careful control of cuticle biosynthesis could produce drought stress tolerance in relevant crop plants[2]. In fleshy fruits, the cuticle maintains adequate water content during fruit development on the plant and reduces water loss in fruit during postharvest[3]. Efforts to elucidate the molecular pathway of cuticle biosynthesis have been carried out for fleshy fruits such as tomato (Solanum lycopersicum)[4], apple (Malus domestica)[5], sweet cherry (Prunus avium)[6], mango (Mangifera indica)[7], and pear (Pyrus 'Yuluxiang')[8].

    The plant cuticle is formed by the two main layers cutin and cuticular waxes[3]. Cutin is composed mainly of oxygenated long-chain (LC) fatty acids (FA), which are synthesized by cytochrome p450 (CYP) enzymes. CYP family 86 subfamily A (CYP86A) enzymes carry out the terminal (ω) oxidation of LC-FA[9]. Then, CYP77A carries out the mid-chain oxidation to synthesize the main cutin monomers. In Arabidopsis, AtCYP77A4 and AtCYP77A6 carry out the synthesis of mid-chain epoxy and mid-chain dihydroxy LC-FA, respectively[10,11]. AtCYP77A6 is required for the cutin biosynthesis and the correct formation of floral surfaces[10]. The expression of CYP77A19 (KF410855) and CYP77A20 (KF410856) from potato (Solanum tuberosum) restored the petal cuticular impermeability in Arabidopsis null mutant cyp77a6-1, tentatively by the synthesis of cutin monomers[12]. In eggplant (Solanum torvum), the over-expression of StoCYP77A2 leads to resistance to Verticillium dahlia infection in tobacco plants[13]. Although the function of CYP77A2 in cutin biosynthesis has not yet been tested, gene expression analysis suggests that CaCYP77A2 (A0A1U8GYB0) could play a role in cutin biosynthesis during pepper fruit development[14].

    It has been hypothesized that the export of cuticle precursors is carried out by ATP binding cassette subfamily G (ABCG) transporters. ABCG11/WBC11, ABCG12, and ABCG13 are required for the load of cuticle lipids in Arabidopsis[15−17], but ABCG13 function appears to be specific to the flower epidermis[18]. The overexpression of TsABCG11 (JQ389853) from Thellungiella salsugineum increases cuticle amounts and promotes tolerance to different abiotic stresses in Arabidopsis[19].

    Once exported, the cutin monomers are polymerized on the surface of epidermal cells. CD1 code for a Gly-Asp-Ser-Leu motif lipase/esterase (GDSL) from tomato required for the cutin formation through 2-mono(10,16-dihydroxyhexadecanoyl)glycerol esterification[20]. GDSL1 from tomato carries out the ester bond cross-links of cutin monomers located at the cuticle layers and is required for cuticle deposition in tomato fruits[21]. It has been shown that the transcription factor MIXTA-like reduces water loss in tomato fruits through the positive regulation of the expression of CYP77A2, ABCG11, and GDSL1[22]. Despite the relevant role of cuticles in maintaining cactus homeostasis in desert environments[1], the molecular mechanism of cuticle biosynthesis has yet to be described for cactus fruits.

    Stenocereus thurberi is a columnar cactus endemic from the Sonoran desert (Mexico), which produces an ovoid-globose fleshy fruit named sweet pitaya[23]. In its mature state, the pulp of sweet pitaya contains around 86% water with a high content of antioxidants and natural pigments such as betalains and phenolic compounds, which have nutraceutical and industrial relevance[23]. Due to the arid environment in which pitaya fruit grows, studying its molecular mechanism of cuticle biosynthesis can generate new insights into understanding species' adaptation mechanisms to arid environments. Nevertheless, sequences of transcripts from S. thurberi in public databases are scarce.

    RNA-sequencing technology (RNA-seq) allows the massive generation of almost all the transcripts from non-model plants, even if no complete assembled genome is available[24]. Recent advances in bioinformatic tools has improved our capacity to identify long non-coding RNA (lncRNA), which have been showed to play regulatory roles in relevant biological processes, such as the regulation of drought stress tolerance in plants[25], fruit development, and ripening[26−29].

    In this study, RNA-seq data were obtained for the de novo assembly and characterization of the S. thurberi fruit peel transcriptome. As a first approach, three transcripts, StCYP77A, StABCG11, and StGDSL1, tentatively involved in cuticle biosynthesis, were identified and quantified during sweet pitaya fruit development. Due to no gene expression analysis having been carried out yet for S. thurberi, stably expressed constitutive genes were identified for the first time.

    Materials and methods

    Plant materials and gene library sequencing

    Sweet pitaya fruits (S. thurberi) without physical damage were hand harvested from plants in a native conditions field located at Carbó, Sonora, México. They were collocated in a cooler containing dry ice and transported immediately to the laboratory. The superficial part of the peels (~1 mm deep) was removed carefully from the fruits using a scalpel. Peel samples from three fruits were pooled according to their tentative stage of development defined by their visual characteristics, frozen in liquid nitrogen, and pulverized to create a single biological replicate. Four samples belonging to four different plants were analyzed. All fruits harvested were close to the ripening stage. Samples named M1 and M2 were turning from green to ripe [~35−40 Days After Flowering (DAF)], whereas samples M3 and M4 were turning from ripe to overripe (~40−45 DAF).

    Total RNA was isolated from the peels through the Hot Borate method[30]. The concentration and purity of RNA were determined in a spectrophotometer Nanodrop 2000 (Thermo Fisher) by measuring the 260/280 and 260/230 absorbance ratios. RNA integrity was evaluated through electrophoresis in agarose gel 1% and a Bioanalyzer 2100 (Agilent). Pure RNA was sequenced in the paired-end mode in an Illumina NextSeq 500 platform at the University of Arizona Genetics Core Facility. Four RNA-seq libraries, each of them from each sample, were obtained, which include a total of 288,199,704 short reads with a length of 150 base pairs (bp). The resulting sequence data can be accessed at the Sequence Read Archive (SRA) repository of the NCBI through the BioProject ID PRJNA1030439. Libraries are named corresponding to the names of samples M1, M2, M3, and M4.

    de novo transcriptome assembly and quality analysis

    FastQC software (www.bioinformatics.babraham.ac.uk/projects/fastqc) was used for short reads quality analysis. Short reads with poor quality were trimmed or eliminated by Trimmomatic (www.usadellab.org/cms/?page=trimmomatic) with a trailing and leading of 25, a sliding window of 4:25, and a minimum read length of 80 bp. A total of 243,194,888 reads with at least a 25 quality score on the Phred scale were used to carry out the de novo assembly by Trinity (https://github.com/trinityrnaseq/trinityrnaseq/wiki) with the following parameters: minimal k-mer coverage of 1, normalization of 50, and minimal transcript length of 200 bp.

    Removal of contaminating sequences and ribosomal RNA (rRNA) was carried out through SeqClean. To remove redundancy, transcripts with equal or more than 90% of identity were merged through CD-hit (www.bioinformatics.org/cd-hit/). Alignment and quantification in terms of transcripts per million (TPM) were carried out through Bowtie (https://bowtie-bio.sourceforge.net/index.shtml) and RSEM (https://github.com/deweylab/RSEM), respectively. Transcripts showing a low expression (TPM < 0.01) were discarded. Assembly quality was evaluated by calculating the parameters N50 value, mean transcript length, TransRate score, and completeness. The statistics of the transcriptome were determined by TrinityStats and TransRate (https://hibberdlab.com/transrate/). The transcriptome completeness was determined through a BLASTn alignment (E value < 1 × 10−3) by BUSCO (https://busco.ezlab.org/) against the database of conserved orthologous genes from Embryophyte.

    Functional annotation of protein-coding transcripts

    To predict the proteins tentatively coded in the S. thurberi transcriptome, the best homology match of the assembled transcripts was found by alignment to the Swiss-Prot, RefSeq, nr-NCBI, PlantTFDB, iTAK, TAIR, and ITAG databases using the BLAST algorithm with an E value threshold of 1 × 10−10 for the nr-NCBI database and of 1 × 10−5 for the others[31−34]. An additional alignment was carried out to the protein databases of commercial fruits Persea americana, Prunus persica, Fragaria vesca, Citrus cinensis, and Vitis vinifera to proteins of the cactus Opuntia streptacantha, and the transcriptomes of the cactus Hylocereus polyrhizus, Pachycereus pringlei, and Selenicereus undatus. The list of all databases and the database websites of commercial fruits and cactus are provided in Supplementary Tables S1 & S2. The open reading frame (ORF) of the transcripts and the protein sequences tentative coded from the sweet pitaya transcriptome was predicted by TransDecoder (https://github.com/TransDecoder/TransDecoder/wiki), considering a minimal ORF length of 75 amino acids (aa). The search for protein domains was carried out by the InterPro database (www.ebi.ac.uk/interpro). Functional categorization was carried out by Blast2GO based on GO terms and KEGG metabolic pathways[35].

    Identification of long non-coding transcripts

    LncRNA were identified based on the methods reported in previous studies[25,29,36]. Transcripts without homology to any protein from Swiss-Prot, RefSeq, nr-NCBI, PlantTFDB, iTAK, TAIR, ITAG, P. americana, P. persica, F. vesca, C. cinensis, V. vinifera, and O. streptacantha databases, without a predicted ORF longer than 75 aa, and without protein domains in the InterPro database were selected to identify tentative lncRNA.

    Transcripts coding for signal peptide or transmembrane helices were identified by SignalP (https://services.healthtech.dtu.dk/services/SignalP-6.0/) and TMHMM (https://services.healthtech.dtu.dk/services/TMHMM-2.0/), respectively, and discarded. Further, transcripts corresponding to other non-coding RNAs (ribosomal RNA and transfer RNA) were identified through Infernal by using the Rfam database[37] and discarded. The remaining transcripts were analyzed by CPC[38], and CPC2[39] to calculate their coding potential. Transcripts with a coding potential score lower than −1 for CPC and a coding probability lower than 0.1 for CPC2 were considered lncRNA. To characterize the identified lncRNA, the length and abundance of coding and lncRNA were calculated. Bowtie and RSEM were used to align and quantify raw counts, respectively. The edgeR package[40] was used to normalize raw count data in terms of counts per million (CPM) for both coding and lncRNA.

    Identification of tentative reference genes

    To obtain the transcript's expression, the aligning of short reads and quantifying of transcripts were carried out through Bowtie and RSEM software, respectively. A differential expression analysis was carried out between the four libraries by edgeR package in R Studio. Only the transcripts with a count equal to or higher than 0.5 in at least one sample were retained for the analysis. Transcripts with log2 Fold Change (log2FC) between +1 and −1 and with a False Discovery Rate (FDR) lower than 0.05 were taken as not differentially expressed (NDE).

    For the identification of the tentative reference genes two strategies were carried out as described below: i) The NDE transcripts were aligned by BLASTn (E value < 1 × 10−5) to 43 constitutive genes previously reported in fruits from the cactus H. polyrhizus, S. monacanthus, and S. undatus[41−43] to identify possible homologous constitutive genes in S. thurberi. Then, the homologous transcripts with the minimal coefficient of variation (CV) were selected; ii) For all the NDE transcripts, the percentile 95 value of the mean CPM and the percentile 5 value of the CV were used as filters to recover the most stably expressed transcripts, based on previous studies[44]. Finally, transcripts to be tested by quantitative reverse transcription polymerase chain reaction (qRT-PCR) were selected based on their homology and tentative biological function.

    Evaluation of reliable reference genes

    The fruit harvesting was carried out as described above. Sweet pitaya fruit takes about 43 d to ripen, therefore, open flowers were tagged, and fruits with 10, 20, 30, 35, and 40 DAF were collected to cover the pitaya fruit development process (Supplementary Fig. S1). The superficial part of the peels (~1 mm deep) was removed carefully from the fruits using a scalpel. Peel samples from three fruits were pooled according to their stage of development defined by their DAF, frozen in liquid nitrogen, and pulverized to create a single biological replicate. One biological replicate consisted of peels from three fruits belonging to the same plant. Two to three biological replicates were evaluated for each developmental stage. Two technical replicates were analyzed for each biological replicate. RNA extraction, quantification, RNA purity, and RNA integrity analysis were carried out as described above.

    cDNA was synthesized from 100 ng of RNA by QuantiTect Reverse Transcription Kit (QIAGEN). Primers were designed using the PrimerQuest™, UNAFold, and OligoAnalyzer™ tools from Integrated DNA Technologies (www.idtdna.com/pages) and following the method proposed by Thornton & Basu[45]. Transcripts quantification was carried out in a QIAquant 96 5 plex according to the PowerUp™ SYBR™ Green Master Mix protocol (Applied Biosystems), with a first denaturation step for 2 min at 95 °C, followed by 40 cycles of denaturation step at 95 °C for 15 s, annealing and extension steps for 30 s at 60 °C.

    The Cycle threshold (Ct) values obtained from the qRT-PCR were analyzed through the algorithms BestKeeper, geNorm, NormFinder, and the delta Ct method[46]. RefFinder (www.ciidirsinaloa.com.mx/RefFinder-master/) was used to integrate the stability results and to find the most stable expressed transcripts in sweet pitaya fruit peel during development. The pairwise variation value (Vn/Vn + 1) was calculated through the geNorm algorithm in R Studio software[47].

    Quantification of cuticle biosynthesis-related transcripts

    An alignment of 17 reported cuticle biosynthesis genes from model plants were carried out by BLASTx against the predicted proteins from sweet pitaya. Two additional alignments of 17 charaterized cuticle biosynthesis proteins from model plants against the transcripts and predicted proteins of sweet pitaya were carried out by tBLASTn and BLASTp, respectively. An E value threshold of 1 × 10−5 was used, and the unique best hits were recovered for all three alignments. The sequences of the 17 characterized cuticle biosynthesis genes and proteins from model plants are showed in Supplementary Table S3. The specific parameters and the unique best hits for all the alignments carried out are shown in Supplementary Tables S4−S8.

    Cuticle biosynthesis-related transcripts tentatively coding for a cytochrome p450 family 77 subfamily A (CYP77A), a Gly-Asp-Ser-Leu motif lipase/esterase 1 (GDSL1), and an ATP binding cassette transporter subfamily G member 11 (ABCG11) were identified by best bi-directional hit according to the functional annotation described above. Protein-conserved domains, signal peptide, and transmembrane helix were predicted through InterProScan, SignalP 6.0, and TMHMM, respectively. Alignment of the protein sequences to tentative orthologous of other plant species was carried out by the MUSCLE algorithm[48]. A neighbor-joining (NJ) phylogenetic tree with a bootstrap of 1,000 replications was constructed by MEGA11[49].

    Fruit sampling, primer design, RNA extraction, cDNA synthesis, and transcript quantification were performed as described above. Relative expression was calculated according to the 2−ΔΔCᴛ method[50]. The sample corresponding to 10 DAF was used as the calibrator. The transcripts StEF1a, StTUA, StUBQ3, and StEF1a + StTUA were used as normalizer genes.

    Statistical analysis

    Normality was assessed according to the Shapiro-Wilk test. Significant differences in the expression of the cuticle biosynthesis-related transcripts between fruit developmental stages were determined by one-way ANOVA based on a completely randomized sampling design and a Tukey honestly significant difference (HSD) test, considering a p-value < 0.05 as significant. Statistical analysis was carried out through the stats package in R Studio.

    Results

    Short reads and assembly quality

    RNA was extracted from the peels of ripe sweet pitaya fruits (S. thurberi) from plants located in the Sonoran Desert, Mexico. Four cDNA libraries were sequenced in an Illumina NextSeq 500 platform at the University of Arizona Genetics Core Facility. A total of 288,199,704 reads with 150 base pairs (bp) in length were sequenced in paired-end mode. After trimming, 243,194,888 (84.38%) cleaned short reads with at least 29 mean quality scores per read in the Phred scale and between 80 to 150 bp in length were obtained to carry out the assembly. After removing contaminating sequences, redundancy, and low-expressed transcripts, the assembly included 174,449 transcripts with an N50 value of 2,110 bp. Table 1 shows the different quality variables of the S. thurberi fruit peel transcriptome. BUSCO score showed that 85.4% are completed transcripts, although out of these, 37.2% were found to be duplicated. The resulting sequence data can be accessed at the SRA repository of the NCBI through the BioProject ID PRJNA1030439.

    Table 1.  Quality metrics of the Stenocereus thurberi fruit peel transcriptome.
    Metric Data
    Total transcripts 174,449
    N50 2,110
    Smallest transcript length (bp) 200
    Largest transcript length (bp) 19,114
    Mean transcript length (bp) 1,198.69
    GC (%) 41.33
    Total assembled bases 209,110,524
    TransRate score 0.05
    BUSCO score (%) C: 85.38 (S:48.22, D:37.16),
    F: 10.69, M: 3.93.
    Values were calculated through the TrinityStats function of Trinity and TransRate software. Completeness analysis was carried out through BUSCO by aligning the transcriptome to the Embryophyte database through BLAST with an E value threshold of 1 × 10−3. Complete (C), single (S), duplicated (D), fragmented (F), missing (M).
     | Show Table
    DownLoad: CSV

    Homology searches

    A summary of the homology search in the main public protein database for the S. thurberi transcriptome is shown in Supplementary Table S1. From these databases, the higher homologous transcripts were found in RefSeq with 93,993 (53.87 %). Based on the E value distribution, for 41,685 (44%) and 68,853 (49%) of the hits, it was found a strong homology (E value lower than 1 × 10−50) to proteins in the Swiss-Prot and RefSeq databases, respectively (Supplementary Fig. S2a & b). On the other hand, 56,539 (52.34%) and 99,599 (71.11%) of the matches showed a percentage of identity higher than 60% in the Swiss-Prot and RefSeq databases, respectively (Supplementary Fig. S2c & d).

    Figure 1 shows the homology between transcripts from S. thurberi and proteins of commercial fruits, as well as proteins and transcripts of cacti. Transcripts from S. thurberi homologous to proteins from fruits of commercial interest avocado (P. americana), peach (P. persica), strawberry (F. vesca), orange (C. sinensis), and grapefruit (V. vinifera) ranged from 77,285 (44.30%) to 85,421 (48.96%), with 70,802 transcripts homologous to all the five fruit protein databases (Fig. 1a).

    Transcripts homologous to transcripts or proteins from the cactus dragon fruit (H. polyrhizus), prickly pear cactus (O. streptacantha), Mexican giant cardon (P. pringlei), and pitahaya (S. undatus) ranged from 76,238 (43.70%) to 114,933 (65.88%), with 64,009 transcripts homologous to all the four cactus databases (Fig. 1b). Further, out of the total of transcripts, 44,040 transcripts (25.25%) showed homology only to sequences from cactus, but not for model plants Arabidopsis, tomato, or the commercial fruits included in this study (Fig. 1c).

    Figure 1.  Venn diagram of the homology search results against model plants databases, commercial fruits, and cactus. The number in the diagram corresponds to the number of transcripts from S. thurberi homologous to sequences from that plant species. (a) Homologous to sequences from Fragaria vesca (Fa), Persea americana (Pa), Prunus persica (Pp), Vitis vinifera (Vv), and Citrus sinensis (Cs). (b) Homologous to sequences from Opuntia streptacantha (Of), Selenicereus undatus (Su), Hylocereus polyrhizus (Hp), and Pachycereus pringlei (Pap). (c) Homologous to sequences from Solanum lycopersicum (Sl), Arabidopsis thaliana (At), from the commercial fruits (Fa, Pa, Pp, Vv, and Cs), or the cactus included in this study (Of, Su, Hp, and Pap). Homologous searching was carried out by BLAST alignment (E value < 1 × 10−5). The Venn diagrams were drawn by ggVennDiagram in R Studio.
    DownLoad: Full-Size Img PowerPoint

    A total of 45,970 (26.35%), 58,704 (33.65%), and 48,186 (27.65%) transcripts showed homology to transcription factors, transcriptional regulators, and protein kinases in the PlantTFDB, iTAK-TR, and iTAK-PK databases, respectively (Supplementary Tables S1, S9−S11). For the PlantTFDB, the homologous transcripts belong to 57 transcriptional factors (TF) families (Fig. 2 & Supplementary Table S9), from which, the most frequent were the basic-helix-loop-helix (bHLH), myeloblastosis-related (MYB-related), NAM, ATAF, and CUC (NAC), ethylene responsive factor (ERF), and the WRKY domain families (WRKY) (Fig. 2).

    Figure 2.  Transcription factor (TF) families distribution of S. thurberi fruit peel transcriptome. The X-axis indicates the number of transcripts with hits to each TF family. Alignment to the PlantTFDB database by BLASTx was carried out with an E value threshold of 1 × 10−5. The bar graph was drawn by ggplot2 in R Studio.
    DownLoad: Full-Size Img PowerPoint

    Functional categorization

    Based on the homology found and the functional domain searches, gene ontology terms (GO) were assigned to 68,559 transcripts (Supplementary Table S12). Figure 3 shows the top 20 GO terms assigned to the S. thurberi transcriptome, corresponding to the Biological Processes (BP) and Molecular Function (MF) categories. For BP, organic substance metabolic processes, primary metabolic processes, and cellular metabolic processes showed a higher number of transcripts (Supplementary Table S13). Further, for MF, organic cyclic compound binding, heterocyclic compound binding, and ion binding were the processes with the higher number of transcripts. S. thurberi transcripts were classified into 142 metabolic pathways from the KEGG database (Supplementary Table S14). The pathways with the higher number of transcripts recorded were pyruvate metabolism, glycerophospholipid metabolism, glycolysis/gluconeogenesis, and citrate cycle. Further, among the top 20 KEEG pathways, the cutin, suberin, and wax biosynthesis include more than 30 transcripts (Fig. 4).

    Figure 3.  Top 20 Gene Ontology (GO) terms assigned to the S. thurberi fruit peel transcriptome. Bars indicate the number of transcripts assigned to each GO term. Assignment of GO terms was carried out by Blast2GO with default parameters. BP and MF mean Biological Processes and Molecular Functions GO categories, respectively. The graph was drawn by ggplot2 in R Studio.
    DownLoad: Full-Size Img PowerPoint
    Figure 4.  Top 20 KEGG metabolic pathways distribution in the S. thurberi fruit peel transcriptome. Bars indicate the number of transcripts assigned to each KEGG pathway. Assignment of KEGG pathways was carried out in the Blast2GO suite. The bar graph was drawn by ggplot2 in R Studio.
    DownLoad: Full-Size Img PowerPoint

    Identification of lncRNA

    Out of the total of transcripts, 43,391 (24.87%) were classified as lncRNA (Supplementary Tables S15 & S16). Figure 5 shows a comparison of the length (Fig. 5a) and expression (Fig. 5b) of lncRNA and coding RNA. Both length and expression values were higher in coding RNA than in lncRNA. In general, coding RNA ranged from 201 to 18,629 bp with a mean length of 1,507.18, whereas lncRNA ranged from 200 to 5,198 bp with a mean length of 481.51 (Fig. 5a). The higher expression values recorded from coding RNA and lncRNA were 12.83 and 9.45 log2(CPM), respectively (Fig. 5b).

    Figure 5.  Comparison of coding RNA and long non-coding RNA (lncRNA) from S. thurberi transcriptome. (a) Box plot of transcript length distribution. The Y-axis indicates the length of each transcript in base pairs. (b) Box plot of expression levels. The Y-axis indicates the log2 of the count per million of reads (log2(CPM)) recorded for each transcript. Expression levels were calculated by the edgeR package in R studio. (a), (b) The lines inside the boxes indicate the median. The higher and lower box limits represent the 75th and 25th percentiles, respectively. The box plots were drawn by ggplot2 in R Studio.
    DownLoad: Full-Size Img PowerPoint

    Identification of tentative reference genes

    To identify the transcripts without significant changes in expression between the four RNA-seq libraries, a differential expression analysis was carried out. Of the total of transcripts, 4,980 were not differentially expressed (NDE) at least in one paired comparison between the libraries (Supplementary Tables S17−S20). Mean counts per million of reads (CPM) and coefficient of variation (CV)[44] were calculated for these NDE transcripts. Transcripts with a CV value lower than 0.113, corresponding with the percentile 5 of the CV, and a mean CPM higher than 1,138.06, corresponding with the percentile 95 of the mean CPM were used as filters to identify the most stably expressed transcripts (Supplementary Table S21). Based on its homology and its tentative biological function, five transcripts were selected to be tested as tentative reference genes. Besides, three NDE transcripts homologous to previously identified stable expressed reference genes in other species of cactus fruit[41−43] were selected (Supplementary Table S22). Homology metrics for the eight tentative reference genes selected are shown in Supplementary Table S23. The primer sequences used to amplify the transcripts by qRT-PCR and their nucleotide sequence are shown in Supplementary Tables S24 & S25, respectively.

    Expression stability of tentative reference genes

    The amplification specificity of the eight candidate reference genes determined by melting curves analysis is shown in Supplementary Fig. S3. For the eight tentative reference transcripts selected, the cycle threshold (Ct) values were recorded during sweet pitaya fruit development by qRT-PCR (Supplementary Table S26). The Ct values obtained ranged from 16.85 to 30.26 (Fig. 6a). Plastidic ATP/ADP-transporter (StTLC1) showed the highest Ct values with a mean of 27.34 (Supplementary Table S26). Polyubiquitin 3 (StUBQ3) showed the lowest Ct values in all five sweet pitaya fruit developmental stages (Fig. 6a).

    Figure 6.  Expression stability analysis of tentative reference genes. (a) Box plot of cycle threshold (Ct) distribution of candidate reference genes during sweet pitaya fruit development (10, 20, 30, 35, and 40 d after flowering). The black line inside the box indicates the median. The higher and lower box limits represent the 75th and 25th percentiles, respectively. (b) Bar chart of the geometric mean (geomean) of ranking values calculated by RefFinder for each tentative reference gene (X-axis). The lowest values indicate the best reference genes. (c) Bar chart of the pairwise variation analysis and determination of the optimal number of reference genes by the geNorm algorithm. A pairwise variation value lower than 0.15 indicates that the use of Vn/Vn + 1 reference genes is reliable for the accurate normalization of qRT-PCR data. The Ct data used in the analysis were calculated by qRT-PCR in a QIAquant 96 5 plex (QIAGEN) according to the manufacturer's protocol. The box plot and the bar graphs were drawn by ggplot2 and Excel programs, respectively. Abbreviations: Actin 7 (StACT7), alpha-tubulin (StTUA), elongation factor 1-alpha (StEF1a), COP1-interactive protein 1 (StCIP1), plasma membrane ATPase 4 (StPMA4), BEL1-like homeodomain protein 1 (StBLH1), polyubiquitin 3 (StUBQ3), and plastidic ATP/ADP-transporter (StTLC1).
    DownLoad: Full-Size Img PowerPoint

    The best stability values calculated by NormFinder were 0.45, 0.51, 0.97, and 0.99, corresponding to the transcripts elongation factor 1-alpha (StEF1a), alpha-tubulin (StTUA), plastidic ATP/ADP-transporter (StTLC1), and actin 7 (StACT7), respectively (Supplementary Table S27). For BestKeeper, the most stable expressed transcripts were StUBQ3, StTUA, and StEF1a, with values of 0.72, 0.75, and 0.87, respectively. In the case of the delta Ct method[51], the transcripts StEF1a, StTUA, and StTLC1 showed the best stability.

    According to geNorm analysis, the most stable expressed transcripts were StTUA, StEF1a, StUBQ3, and StACT7, with values of 0.74, 0.74, 0.82, and 0.96, respectively. All the pairwise variation values (Vn/Vn + 1) were lower than 0.15, ranging from 0.019 for V2/V3 to 0.01 for V6/V7 (Fig. 6c). The V value of 0.019 obtained for V2/V3 indicates that the use of the best two reference genes StTUA and StEF1a is reliable enough for the accurate normalization of qRT-PCR data, therefore no third reference gene is required[47]. Except for BestKeeper analysis, StEF1a and StTUA were the most stable transcripts for all of the methods carried out in this study (Supplementary Table S27). The comprehensive ranking analysis indicates that StEF1a, followed by StTUA and StUBQ3, are the most stable expressed genes and are stable enough to be used as reference genes in qRT-PCR analysis during sweet pitaya fruit development (Fig. 6b).

    Identification of cuticle biosynthesis-related transcripts

    Three cuticle biosynthesis-related transcripts TRINITY_DN17030_c0_g1_i2, TRINITY_DN15394_c0_g1_i1, and TRINITY_DN23528_c1_g1_i1 tentatively coding for the enzymes cytochrome p450 family 77 subfamily A (CYP77A), Gly-Asp-Ser-Leu motif lipase/esterase 1 (GDSL1), and an ATP binding cassette transporter subfamily G member 11 (ABCG11/WBC11), respectively, were identified and quantified. The nucleotide sequence and predicted amino acid sequences of the three transcripts are shown in Supplementary File 1. The best homology match for StCYP77A (TRINITY_DN17030_c0_g1_i2) was for AtCYP77A4 (AT5G04660) from Arabidopsis and SmCYP77A2 (P37124) from eggplant (Solanum melongena) in the TAIR and Swiss-Prot databases, respectively (Supplementary Table S23).

    TransDecoder, InterPro, and TMHMM analysis showed that StCYP77A codes a polypeptide of 518 amino acids (aa) in length that comprises a cytochrome P450 E-class domain (IPR002401) and a transmembrane region (residues 10 to 32). The phylogenetic tree constructed showed that StCYP77A is grouped in a cluster with all the CYP77A2 proteins included in this analysis, being closer to CYP77A2 (XP_010694692) from B. vulgaris and Cgig2_012892 (KAJ8441854) from Carnegiea gigantean (Supplementary Fig. S4).

    StGDSL1 (TRINITY_DN15394_c0_g1_i1) alignment showed that it is homologous to a GDSL esterase/lipase from Arabidopsis (Q9LU14) and tomato (Solyc03g121180) (Supplementary Table S23). TransDecoder, InterPro, and SignalP analysis showed that StGDSL1 codes a polypeptide of 354 aa in length that comprises a GDSL lipase/esterase domain IPR001087 and a signal peptide with a cleavage site between position 25 and 26 (Supplementary Fig. S5).

    Supplementary Figure S6 shows the analysis carried out on the predicted amino acid sequence of StABCG11 (TRINITY_DN23528_c1_g1_i1). The phylogenetic tree constructed shows three clades corresponding to the ABCG13, ABCG12, and ABCG11 protein classes with bootstrap support ranging from 40% to 100% (Supplementary Fig. S6a). StABCG11 is grouped with all the ABCG11 transporters included in this study in a well-separated clade, being closely related to its tentative ortholog from C. gigantean Cgig2_004465 (KAJ8441854). InterPro and TMHMM results showed that the StABCG11 sequence contains an ABC-2 type transporter transmembrane domain (IPR013525; PF01061.27) with six transmembrane helices (Supplementary Fig. S6b).

    The predicted protein sequence of StABCG11 is 710 aa in length, holding the ATP binding domain (IPR003439; PF00005.30) and the P-loop containing nucleoside triphosphate hydrolase domain (IPR043926; PF19055.3) of the ABC transporters of the G family. Multiple sequence alignment shows that the Walker A and B motif sequence and the ABC signature[15] are conserved between the ABCG11 transporters from Arabidopsis, tomato, S. thurberi, and C. gigantean (Supplementary Fig. S6c).

    Evaluation of reliable reference genes and quantification of cuticle biosynthesis-related transcripts

    According to the results of the expression stability analysis (Fig. 6), four normalization strategies were tested to quantify the three cuticle biosynthesis-related transcripts during sweet pitaya fruit development. The four strategies consist of normalizing by StEF1a, StTUA, StUBQ3, or StEF1a+StTUA. Primer sequences used to quantify the transcripts StCYP77A (TRINITY_DN17030_c0_g1_i2), StGDSL1 (TRINITY_DN15394_c0_g1_i1), and StABCG11 (TRINITY_DN23528_c1_g1_i1) by qRT-PCR during sweet pitaya fruit development are shown in Supplementary Table S24.

    The three cuticle biosynthesis-related transcripts showed differences in expression during sweet pitaya fruit development (Supplementary Table S28). The same expression pattern was recorded for the three cuticle biosynthesis transcripts when normalization was carried out by StEF1a, StTUA, StUBQ3, or StEF1a + StTUA (Fig. 7). A higher expression of StCYP77A and StGDSL1 are shown at the 10 and 20 DAF, showing a decrease at 30, 35, and 40 DAF. StABCG11 showed a similar behavior, with a higher expression at 10 and 20 DAF and a reduction at 30 and 35 DAF. Nevertheless, unlike StCYP77A and StGDSL1, a significant increase at 40 DAF, reaching the same expression as compared with 10 DAF, is shown for StABCG11 (Fig. 7).

    Figure 7.  Expression analysis of cuticle biosynthesis-related transcripts StCYP77A, StGDSL1, and StABCG11 during sweet pitaya (Stenocereus thurberi) fruit development. Relative expression was calculated through the 2−ΔΔCᴛ method using elongation factor 1-alpha (StEF1a), alpha-tubulin (StTUA), polyubiquitin 3 (StUBQ3), or StEF1a + StTUA as normalizing genes at 10, 20, 30, 35, and 40 d after flowering (DAF). The Y-axis and error bars represent the mean of the relative expression ± standard error (n = 4−6) for each developmental stage in DAF. The Ct data for the analysis was recorded by qRT-PCR in a QIAquant 96 5 plex (QIAGEN) according to the manufacturer's protocol. The graph line was drawn by ggplot2 in R Studio. Abbreviations: cytochrome p450 family 77 subfamily A (StCYP77A), Gly-Asp-Ser-Leu motif lipase/esterase 1 (StGDSL1), and ATP binding cassette transporter subfamily G member 11 (StABCG11).
    DownLoad: Full-Size Img PowerPoint

    Discussion

    S. thurberi transcriptome quality is similar to that reported for de novo assembly of non-model plants

    Characteristics of a well-assembled transcriptome include an N50 value closer to 2,000 bp, a high percentage of conserved transcripts completely assembled (> 80%), and a high proportion of reads mapping back to the assembled transcripts[52]. In the present study, the first collection of 174,449 transcripts from S. thurberi fruit peel are reported. The generated transcriptome showed an N50 value of 2,110 bp, a TransRate score of 0.05, and a GC percentage of 41.33 (Table 1), similar to that reported for other de novo plant transcriptome assemblies[53]. According to BUSCO, 85.4% of the orthologous genes from the Embryophyta databases completely matched the S. thurberi transcriptome, and only 3.9% were missing (Table 1). These results show that the S. thurberi transcriptome generated is not fragmented, and it is helpful in predicting the sequence of almost all the transcripts expressed in sweet pitaya fruit peel[24].

    S. thurberi transcriptome shows higher homology to sequences from species of the same family

    The percentage of transcripts homologous found, E values, and identity distribution (Supplementary Tables S1 & S2; Supplementary Fig. S2) were similar to that reported in the de novo transcriptome assembly for non-model plants and other cactus fruits[41−43,54] and further suggests that the transcriptome assembled of S. thurberi peel is robust[52]. Of the total of transcripts, 70,802 were common to all the five commercial fruit protein databases included in this study, which is helpful for the search for conserved orthologous involved in fruit development and ripening (Fig. 2a). A total of 34,513 transcripts (20%) show homology only to sequences in the cactus's databases, but not in the others (Supplementary Tables S1 & S2; Fig. 1c). This could suggest that a significant conservation of sequences among plants of the Cactaceae family exists which most likely are to have a function in this species adaptation to desert ecosystems.

    Functional categorization analysis suggests a response to stress and an active cuticle biosynthesis in fruit pitaya peel

    To infer the biological functionality represented by the S. thurberi fruit peel transcriptome, gene ontology (GO) terms and KEGG pathways were assigned. Of the main metabolic pathways assigned, 'glycerolipid metabolism' and 'cutin, suberine, and wax biosynthesis' suggests an active cuticle biosynthesis in pitaya fruit peel (Fig. 4). In agreement with the above, the main GO terms assigned for the molecular function (MF) category were 'organic cyclic compound binding', 'transmembrane transporter activity', and 'lipid binding' (Fig. 3). For the biological processes (BP) category, the critical GO terms for the present research are 'cellular response to stimulus', 'response to stress', 'anatomical structure development', and 'transmembrane transport', which could suggest the active development of the fruit epidermis and cuticle biosynthesis for protection to stress.

    The most frequent transcription factors (TF) families found in S. thurberi transcriptome were NAC, WRKY, bHLH, ERF, and MYB-related (Fig. 2), which had been reported to play a function in the tolerance to abiotic stress in plants[55,56]. Although the role of NAC, WRKY, bHLH, ERF, and MYB TF in improving drought tolerance in relevant crop plants has been widely documented[57,58], their contribution to the adaptation of cactus to arid ecosystems has not yet been elucidated and further experimental pieces of evidence are needed.

    It has been reported that the heterologous expression of ERF TF from Medicago truncatula induces drought tolerance and cuticle wax biosynthesis in Arabidopsis leaf[59]. In tomato fruits, the gene SlMIXTA-like which encodes a MYB transcription factor avoids water loss through the positive regulation of genes related to the biosynthesis and transport of cuticle compounds[22]. Despite the relevant role of cuticles in maintaining cactus physiology in desert environments, experimental evidence showing the role of the different TF-inducing cuticle biosynthesis has yet to be reported for cactus fruits.

    lncRNA from pitaya are similar in length and expression to that reported from other plants

    Out of the transcripts, 43,391 were classified as lncRNA (Supplementary Tables S15 & S16). This is the first report of lncRNA identification for the species S. thurberi. In fruits, 3,679 lncRNA has been identified from tomato[26], 3,330 from peach (P. persica)[29], 3,857 from melon (Cucumis melo)[28], 2,505 from hot pepper (Capsicum annuum)[27], and 3,194 from pomegranate (Punica granatum)[36]. Despite the stringent criteria to classify the lncRNA of sweet pitaya fruit (S. thurberi), a higher number of lncRNAs are shown when compared with previous reports. This finding is most likely due to the higher level of redundancy found during the transcriptome analysis. To reduce this redundancy, further efforts to achieve the complete genome assembly of S. thurberi are needed.

    Previous studies showed that lncRNA is shorter and has lower expression levels than coding RNA[60−62]. In agreement with those findings, both the length and expression values of lncRNA from S. thurberi were lower than coding RNA (Fig. 5). It has been suggested that lncRNA could be involved in the biosynthesis of cuticle components in cabbage[61] and pomegranate[36] and that they could be involved in the tolerance to water deficit through the regulation of cuticle biosynthesis in wild banana[60]. Nevertheless, the molecular mechanism by which lncRNA may regulate the cuticle biosynthesis in S. thurberi fruits has not yet been elucidated.

    Previously reported housekeeping genes were found to have the most stable expression during pitaya fruit development

    A relatively constant level of expression characterizes housekeeping genes because they are involved in essential cellular functions. These genes are not induced under specific conditions such as biotic or abiotic stress. Because of this, they are very useful as internal reference genes for qRT-PCR data normalization[63]. Nevertheless, their expression could change depending on plant species, developmental stages, and experimental conditions[64]. Reliable reference genes for a specific experiment in a given species must be identified to carry out an accurate qRT-PCR data normalization[63]. An initial screening of the transcript expression pattern through RNA-seq improves the identification of stably expressed transcripts by qRT-PCR[44,64].

    Identification of stable expressed reference transcripts during fruit development has been carried out in blueberry (Vaccinium bracteatum)[65], kiwifruit (Actinidia chinensis)[66], peach (P. persica)[67], apple (Malus domestica)[68], and soursop (Annona muricata)[69]. These studies include the expression stability analysis through geNorm, NormFinder, and BestKeeper algorithms[68,69], some of which are supported in RNA-seq data[65,66]. Improvement of expression stability analysis by RNA-seq had led to the identification of non-previously reported reference genes with a more stable expression during fruit development than commonly known housekeeping genes in grapevine (V. vinifera)[44], pear (Pyrus pyrifolia and P. calleryana)[64], and pepper (C. annuum)[70].

    For fruits of the Cactaceae family, only a few studies identifying reliable reference genes have been reported[41−43]. Mainly because gene expression analysis has not been carried out previously for sweet pitaya (S. thurberi), the RNA-seq data generated in this work along with geNorm, NormFinder, BestKeeper, and RefFinder algorithms were used to identify reliable reference genes. The comprehensive ranking analysis showed that out of the eight candidate genes tested, StEF1a followed by StTUA and StUBQ3 were the most stable (Fig. 6b). All the pairwise variation values (Vn/Vn + 1) were lower than 0.15 (Fig. 6c), which indicates that StEF1a, StTUA, and StUBQ3 alone or the use of StEF1a and StTUA together are reliable enough to normalize the gene expression data generated by qRT-PCR.

    The genes StEF1a, StTUA, and StUBQ3 are homologous to transcripts found in the cactus species known as dragonfruit (Hylocereus monacanthus and H. undatus)[41], which have been tested as tentative reference genes during fruit development. EF1a has been proposed as a reliable reference gene in the analysis of changes in gene expression of dragon fruit (H. monacanthus and H. undatus)[41], peach (P. persica)[67], apple (M. domestica)[68], and soursop (A. muricata)[69]. According to the expression stability analysis carried out in the present study (Fig. 6) four normalization strategies were designed. The same gene expression pattern was recorded for the three target transcripts evaluated when normalization was carried out by the genes StEF1a, StTUA, StUBQ3, or StEF1a + StTUA (Fig. 7). Further, these data indicates that these reference genes are reliable enough to be used in qRT-PCR experiments during fruit development of S. thurberi.

    Cutin biosynthesis could have a relevant role in the first stages of pitaya fruit development

    The plant cuticle is formed by two main layers: the cutin, composed mainly of mid-chain oxygenated LC fatty acids, and the cuticular wax, composed mainly of very long-chain (VLC) fatty acids, and their derivates VLC alkanes, VLC primary alcohols, VLC ketones, VLC aldehydes, and VLC esters[3]. In Arabidopsis CYP77A4 and CYP77A6 catalyze the synthesis of midchain epoxy and hydroxy ω-OH long-chain fatty acids, respectively[10,11], which are the main components of fleshy fruit cuticle[3].

    The functional domain search carried out in the present study showed that StCYP77A comprises a cytochrome P450 E-class domain (IPR002401) and a membrane-spanning region from residues 10 to 32 (Supplementary Fig. S4). This membrane-spanning region has been previously characterized in CYP77A enzymes from A. thaliana and Brassica napus[11,71]. It suggests that the protein coded by StCYP77A could catalyze the oxidation of fatty acids embedded in the endoplasmic reticulum membrane of the epidermal cells of S. thurberi fruit. Phylogenetic analysis showed that StCYP77A was closer to proteins from its phylogenetic-related species B. vulgaris (BvCYP772; XP_010694692) and C. gigantea (Cgig2_012892) (Supplementary Fig. S4). StCYP77A, BvCYP77A2, and Cgig2_012892 were closer to SlCYP77A2 and SmCYP77A2 than to CYP77A4 and CYP77A6 proteins, suggesting that StCYP77A (TRINITY_DN17030_c0_g1_i2) could correspond to a CYP77A2 protein.

    Five CYP77A are present in the Arabidopsis genome, named CYP77A4, CYP77A5, CYP77A6, CYP77A7, and CYP77A9, but their role in cuticle biosynthesis has only been reported for CYP77A4 and CYP77A6[72]. It has been suggested that CYP77A2 from eggplant (S. torvum) could contribute to the defense against fungal phytopathogen infection by the synthesis of specific compounds[13]. In pepper fruit (C. annuum), the expression pattern of CYP77A2 (A0A1U8GYB0) and ABCG11 (LOC107862760) suggests a role of CYP77A2 and ABCG11 in cutin biosynthesis at the early stages of pepper fruit development[14].

    In the case of the protein encoded by StGDSL1 (354 aa), the length found in this work is similar to the length of its homologous from Arabidopsis (AT3G16370) and tomato (Solyc03g121180) (Supplementary Fig. S5). A GDSL1 protein named CD1 polymerizes midchain oxygenated ω-OH long-chain fatty acids to form the cutin polyester in the extracellular space of tomato fruit peel[20,21]. It has been suggested that the 25-amino acid N-signal peptide found in StGDSL1 (Supplementary Fig. S5), previously reported in GDSL1 from Arabidopsis, B. napus, and tomato, plays a role during the protein exportation to the extracellular space[21,73].

    A higher expression of StCYP77A, StGDSL1, and StABCG11 is shown at the 10 and 20 DAF of sweet pitaya fruit development (Fig. 7), suggesting the active cuticle biosynthesis at the early stages of sweet pitaya fruit development. In agreement with that, two genes coding for GDSL lipase/hydrolases from tomato named SGN-U583101 and SGN-U579520 are highly expressed in the early stages and during the expansion stages of tomato fruit development, but their expression decreases in later stages[74]. It has been shown that the expression of GDSL genes, like CD1 from tomato, is higher in growing fruit[20,21]. Like tomato, the increase in expression of StCYP77A and StGDSL1 shown in pitaya fruit development could be due to an increase in cuticle deposition during the expansion of the fruit epidermis[20].

    StABCG11 could be playing a relevant role in the transport of cuticle components during sweet pitaya fruit ripening

    The phylogenetic analysis, the functional domains, and the six transmembrane helices found in the StABCG11 predicted protein (Supplementary Fig. S6), suggests that it is an ABCG plasma membrane transporter of sweet pitaya[15]. Indeed, an increased expression of StABCG11 at 40 DAF was recorded in the present study (Fig. 7). Further, this data strongly suggests that it could be playing a relevant role in the transport of cuticle components at the beginning and during sweet pitaya fruit ripening.

    In Arabidopsis, ABCG11 (WBC11) exports cuticular wax and cutin compounds from the plasma membrane[15,75]. It has been reported that a high expression of the ABC plasma membrane transporter from mango MiWBC11 correlates with a higher cuticle deposition during fruit development[7]. The expression pattern for StABCG11, StCYP77A, and StGDSL1 suggests a role of StABCG11 as a cutin compound transporter in the earlier stages of sweet pitaya fruit development (Fig. 7). Further, its increase at 40 DAF suggests that it could be transporting cuticle compounds other than oxygenated long-chain fatty acids, or long-chain fatty acids that are not synthesized by StCYP77A and StGDSL1 in the later stages of fruit development.

    Like sweet pitaya, during sweet cherry fruit (Prunus avium) development, the expression of PaWCB11, homologous to AtABCG11 (AT1G17840), increases at the earlier stages of fruit development decreases at the intermediate stages, and increases again at the later stages[76]. PaWCB11 expression correlated with cuticle membrane deposition at the earlier and intermediate stages of sweet cherry fruit development but not at the later[76]. The increased expression of StABCG11 found in the present study could be due to the increased transport of cuticular wax compounds, such as VLC fatty acids and their derivates, in the later stages of sweet pitaya development[15,75].

    Cuticular waxes make up the smallest amount of the fruit cuticle. Even so, they mainly contribute to the impermeability of the fruit's epidermis[3]. An increase in the transport of cuticular waxes at the beginning of the ripening stage carried out by ABCG transporters could be due to a greater need to avoid water loss and to maintain an adequate amount of water during the ripening of the sweet pitaya fruit. Nevertheless, further expression analysis of cuticular wax biosynthesis-related genes, complemented with chemical composition analysis of cuticles could contribute to elucidating the molecular mechanism of cuticle biosynthesis in cacti and their physiological contribution during fruit development.

    Conclusions

    In this study, the transcriptome of the sweet pitaya (S. thurberi) fruit peel was assembled for the first time. The reference genes found here are a helpful tool for further gene expression analysis in sweet pitaya fruit. Transcripts tentatively involved in cuticle compound biosynthesis and transport are reported for the first time in sweet pitaya. The results suggest a relevant role of cuticle compound biosynthesis and transport at the early and later stages of fruit development. The information generated will help to improve the elucidation of the molecular mechanism of cuticle biosynthesis in S. thurberi and other cactus species in the future. Understanding the cuticle's physiological function in the adaptation of the Cactaceae family to harsh environmental conditions could help design strategies to increase the resistance of other species to face the increase in water scarcity for agricultural production predicted for the following years.

    Author contributions

    The authors confirm contribution to the paper as follows: study conception and design: Tiznado-Hernández ME, Tafolla-Arellano JC, García-Coronado H, Hernández-Oñate MÁ; data collection: Tiznado-Hernández ME, Tafolla-Arellano JC, García-Coronado H, Hernández-Oñate MÁ; analysis and interpretation of results: Tiznado-Hernández ME, García-Coronado H, Hernández-Oñate MÁ, Burgara-Estrella AJ; draft manuscript preparation: Tiznado-Hernández ME, García-Coronado H. All authors reviewed the results and approved the final version of the manuscript.

    Data availability

    All data generated or analyzed during this study are included in this published article and its supplementary information files. The sequence data can be accessed at the Sequence Read Archive (SRA) repository of the NCBI through the BioProject ID PRJNA1030439.

    Acknowledgments

    The authors wish to acknowledge the financial support of Consejo Nacional de Humanidades, Ciencias y Tecnologías de México (CONAHCYT) through project number 579: Elucidación del Mecanismo Molecular de Biosíntesis de Cutícula Utilizando como Modelo Frutas Tropicales. We appreciate the University of Arizona Genetics Core and Illumina for providing reagents and equipment for library sequencing. The author, Heriberto García-Coronado (CVU 490952), thanks the CONAHCYT (acronym in Spanish) for the Ph.D. scholarship assigned (749341). The author, Heriberto García-Coronado, thanks Dr. Edmundo Domínguez-Rosas for the technical support in bioinformatics for identifying long non-coding RNA.

  • The authors declare that they have no conflict of interest.

  • Supplemental Table S1 RNA-seq library and assembly statistics.
    Supplemental Table S2 Summary of RNA-seq annotations found in seven databases.
    Supplemental Table S3 Lists of up- and down-regulated DEGs found in ‘Predator’ at 7 and 14 d of heat stress.
    Supplemental Table S4 Lists of up- and down-regulated DEGs found in ‘Reliant IV’ at 7 and 14 d of heat stress.
    Supplemental Table S5 GO term enrichment of DEGs in both ‘Predator’ and ‘Reliant IV’ at 7 and 14 d of heat stress.
    Supplemental Table S6 WGCNA analysis of DEGs showed two modules (MElightcyan and MEdarkmagenta) and pair-wise gene connection weight in each of them.
    Supplemental Table S7 Gene expression pattern of the central hub genes and connecting genes in MElightcyan and MEdarkmagenta modules, showing as log of fold change(logFC) ± standard error(SE, n = 3). N.S. means ‘not significant’.
    Supplemental Fig. S1 Cluster dendrograms of WGCNA analysis for DEGs in ‘Predator’ (a)  and ‘Reliant IV’ (b)  exposed to 7 d and 14 d of heat stress. The resulting original and merged module eigengenes after dynamic tree cut method were also shown at the bottom.
  • [1] Ortiz R, Sayre KD, Govaerts B, Gupta R, Subbarao GV, et al. 2008. Climate change: can wheat beat the heat? Agriculture, Ecosystems & Environment 126:46−58 doi: 10.1016/j.agee.2008.01.019

    CrossRef   Google Scholar

    [2] Change IPOC. 2014. IPCC. Climate change
    [3] Wahid A, Gelani S, Ashraf M, Foolad MR. 2007. Heat tolerance in plants: an overview. Environmental and Experimental Botany 61:199−223 doi: 10.1016/j.envexpbot.2007.05.011

    CrossRef   Google Scholar

    [4] Kotak S, Larkindale J, Lee U, von Koskull-Döring P, Vierling E, et al. 2007. Complexity of the heat stress response in plants. Current Opinion in Plant Biology 10:310−16 doi: 10.1016/j.pbi.2007.04.011

    CrossRef   Google Scholar

    [5] Mittler R, Finka A, Goloubinoff P. 2012. How do plants feel the heat? Trends in Biochemical Sciences 37:118−25 doi: 10.1016/j.tibs.2011.11.007

    CrossRef   Google Scholar

    [6] Qu A, Ding Y, Jiang Q, Zhu C. 2013. Molecular mechanisms of the plant heat stress response. Biochemical and Biophysical Research Communications 432:203−7 doi: 10.1016/j.bbrc.2013.01.104

    CrossRef   Google Scholar

    [7] Bita CE, Gerats T. 2013. Plant tolerance to high temperature in a changing environment: scientific fundamentals and production of heat stress-tolerant crops. Frontiers in Plant Science 4:273 doi: 10.3389/fpls.2013.00273

    CrossRef   Google Scholar

    [8] Ohama N, Sato H, Shinozaki K, Yamaguchi-Shinozaki K. 2017. Transcriptional regulatory network of plant heat stress response. Trends in Plant Science 22:53−65 doi: 10.1016/j.tplants.2016.08.015

    CrossRef   Google Scholar

    [9] Li Q, He Y, Tu M, Yan J, Yu L, et al. 2019. Transcriptome sequencing of two Kentucky bluegrass (Poa pratensis L.) genotypes in response to heat Stress. Notulae Botanicae Horti Agrobotanici Cluj-Napoca 47:328−38 doi: 10.15835/nbha47111365

    CrossRef   Google Scholar

    [10] Wang K, Liu Y, Tian J, Huang K, Shi T, et al. 2017. Transcriptional profiling and identification of heat-responsive genes in perennial ryegrass by RNA-sequencing. Frontiers in Plant Science 8:1032 doi: 10.3389/fpls.2017.01032

    CrossRef   Google Scholar

    [11] Li Z, Cheng B, Zeng W, Liu Z, Peng Y. 2019. The transcriptional and post-transcriptional regulation in perennial creeping bentgrass in response to γ-aminobutyric acid (GABA) and heat stress. Environmental and Experimental Botany 162:515−24 doi: 10.1016/j.envexpbot.2019.03.026

    CrossRef   Google Scholar

    [12] Hu T, Sun X, Zhang X, Nevo E, Fu J. 2014. An RNA sequencing transcriptome analysis of the high-temperature stressed tall fescue reveals novel insights into plant thermotolerance. BMC Genomics 15:1147 doi: 10.1186/1471-2164-15-1147

    CrossRef   Google Scholar

    [13] Li Y, Wang Y, Tang Y, Kakani VG, Mahalingam R. 2013. Transcriptome analysis of heat stress response in switchgrass (Panicum virgatum L.). BMC Plant Biology 13:153 doi: 10.1186/1471-2229-13-153

    CrossRef   Google Scholar

    [14] Xu Y, Wang J, Bonos SA, Meyer WA, Huang B. 2018. Candidate genes and molecular markers correlated to physiological traits for heat tolerance in fine fescue cultivars. International Journal of Molecular Sciences 19:116 doi: 10.3390/ijms19010116

    CrossRef   Google Scholar

    [15] Losvik A, Beste L, Glinwood R, Ivarson E, Stephens J, et al. 2017. Overexpression and down-regulation of barley lipoxygenase LOX2.2 affects jasmonate-regulated genes and aphid fecundity. International Journal of Molecular Sciences 18:2765 doi: 10.3390/ijms18122765

    CrossRef   Google Scholar

    [16] Viswanath KK, Varakumar P, Pamuru RR, Basha SJ, Mehta S, et al. 2020. Plant lipoxygenases and their role in plant physiology. Journal of Plant Biology 63:83−95 doi: 10.1007/s12374-020-09241-x

    CrossRef   Google Scholar

    [17] Umate P. 2011. Genome-wide analysis of lipoxygenase gene family in Arabidopsis and rice. Plant Signaling & Behavior 6:335−38 doi: 10.4161/psb.6.3.13546

    CrossRef   Google Scholar

    [18] Melan MA, Dong X, Endara ME, Davis KR, Ausubel FM, et al. 1993. An Arabidopsis thaliana LIPOXYGENASE1 gene can be induced by pathogens, abscisic acid, and methyl jasmonate. Plant Physiology 101:441−50 doi: 10.1104/pp.101.2.441

    CrossRef   Google Scholar

    [19] Keunen E, Remans T, Opdenakker K, Jozefczak M, Gielen H, et al. 2013. A mutant of the Arabidopsis thaliana LIPOXYGENASE1 gene shows altered signalling and oxidative stress related responses after cadmium exposure. Plant Physiology and Biochemistry 63:272−80 doi: 10.1016/j.plaphy.2012.12.005

    CrossRef   Google Scholar

    [20] Bannenberg G, Martínez M, Hamberg M, Castresana C. 2009. Diversity of the enzymatic activity in the lipoxygenase gene family of Arabidopsis thaliana. Lipids 44:85 doi: 10.1007/s11745-008-3245-7

    CrossRef   Google Scholar

    [21] Ding H, Lai J, Wu Q, Zhang S, Chen L, et al. 2016. Jasmonate complements the function of Arabidopsis lipoxygenase3 in salinity stress response. Plant Science 244:1−7 doi: 10.1016/j.plantsci.2015.11.009

    CrossRef   Google Scholar

    [22] Mack AJ, Peterman TK, Siedow JN. 1987. Lipoxygenase isozymes in higher plants: biochemical properties and physiological role. Isozymes 13:127−54

    Google Scholar

    [23] Hildebrand DF. 1989. Lipoxygenases. Physiologia Plantarum 76:249−53 doi: 10.1111/j.1399-3054.1989.tb05641.x

    CrossRef   Google Scholar

    [24] Siedow JN. 1991. Plant lipoxygenase: structure and function. Annual Review of Plant Physiology and Plant Molecular Biology 42:145−88 doi: 10.1146/annurev.pp.42.060191.001045

    CrossRef   Google Scholar

    [25] Ferrer JL, Austin MB, Stewart C Jr, Noel JP. 2008. Structure and function of enzymes involved in the biosynthesis of phenylpropanoids. Plant Physiology and Biochemistry 46:356−70 doi: 10.1016/j.plaphy.2007.12.009

    CrossRef   Google Scholar

    [26] Huang J, Gu M, Lai Z, Fan B, Shi K, et al. 2010. Functional analysis of the Arabidopsis PAL gene family in plant growth, development, and response to environmental stress. Plant Physiology 153:1526−38 doi: 10.1104/pp.110.157370

    CrossRef   Google Scholar

    [27] Gharibi S, Sayed Tabatabaei BE, Saeidi G, Talebi M, Matkowski A. 2019. The effect of drought stress on polyphenolic compounds and expression of flavonoid biosynthesis related genes in Achillea pachycephala Rech. f. Phytochemistry 162:90−98 doi: 10.1016/j.phytochem.2019.03.004

    CrossRef   Google Scholar

    [28] Chun HJ, Baek D, Cho HM, Lee SH, Jin BJ, et al. 2019. Lignin biosynthesis genes play critical roles in the adaptation of Arabidopsis plants to high-salt stress. Plant Signaling & Behavior 14:1625697 doi: 10.1080/15592324.2019.1625697

    CrossRef   Google Scholar

    [29] Dittrich H, Kutchan TM. 1991. Molecular cloning, expression, and induction of berberine bridge enzyme, an enzyme essential to the formation of benzophenanthridine alkaloids in the response of plants to pathogenic attack. PNAS 88:9969−73 doi: 10.1073/pnas.88.22.9969

    CrossRef   Google Scholar

    [30] Facchini PJ, Penzes C, Johnson AG, Bull D. 1996. Molecular characterization of berberine bridge enzyme genes from opium poppy. Plant Physiology 112:1669−77 doi: 10.1104/pp.112.4.1669

    CrossRef   Google Scholar

    [31] Parani M, Rudrabhatla S, Myers R, Weirich H, Smith B, et al. 2004. Microarray analysis of nitric oxide responsive transcripts in Arabidopsis. Plant Biotechnology Journal 2:359−66 doi: 10.1111/j.1467-7652.2004.00085.x

    CrossRef   Google Scholar

    [32] Locci F, Benedetti M, Pontiggia D, Citterico M, Caprari C, et al. 2019. An Arabidopsis berberine bridge enzyme-like protein specifically oxidizes cellulose oligomers and plays a role in immunity. The Plant Journal 98:540−54 doi: 10.1111/tpj.14237

    CrossRef   Google Scholar

    [33] Santamaría ME, Arnaiz A, Velasco-Arroyo B, Grbic V, Diaz I, et al. 2018. Arabidopsis response to the spider mite Tetranychus urticae depends on the regulation of reactive oxygen species homeostasis. Scientific Reports 8:9432 doi: 10.1038/s41598-018-27904-1

    CrossRef   Google Scholar

    [34] Kiani D, Soltanloo H, Ramezanpour SS, Nasrolahnezhad Qumi AA, Yamchi A, et al. 2017. A barley mutant with improved salt tolerance through ion homeostasis and ROS scavenging under salt stress. Acta Physiologiae Plantarum 39:90 doi: 10.1007/s11738-017-2359-z

    CrossRef   Google Scholar

    [35] Lee RM, Thimmapuram J, Thinglum KA, Gong G, Hernandez AG, et al. 2009. Sampling the waterhemp (Amaranthus tuberculatus) genome using pyrosequencing technology. Weed Science 57:463−69 doi: 10.1614/WS-09-021.1

    CrossRef   Google Scholar

    [36] Gaines TA, Lorentz L, Figge A, Herrmann J, Maiwald F, et al. 2014. RNA-Seq transcriptome analysis to identify genes involved in metabolism-based diclofop resistance in Lolium rigidum. The Plant Journal 78:865−76 doi: 10.1111/tpj.12514

    CrossRef   Google Scholar

    [37] Halkier BA, Møller BL. 1989. Biosynthesis of the cyanogenic glucoside dhurrin in seedlings of Sorghum bicolor (L.) Moench and partial purification of the enzyme system involved. Plant Physiology 90:1552−59 doi: 10.1104/pp.90.4.1552

    CrossRef   Google Scholar

    [38] Sibbesen O, Koch B, Halkier BA, Møller B. 1994. Isolation of the heme-thiolate enzyme cytochrome P-450TYR, which catalyzes the committed step in the biosynthesis of the cyanogenic glucoside dhurrin in Sorghum bicolor (L.) Moench. PNAS 91:9740−44 doi: 10.1073/pnas.91.21.9740

    CrossRef   Google Scholar

    [39] Bak S, Kahn RA, Nielsen HL, Møller BL, Halkier BA. 1998. Cloning of three A-type cytochromes P450, CYP71E1, CYP98, and CYP99 from Sorghum bicolor (L.) Moench by a PCR approach and identification by expression in Escherichia coli of CYP71E1 as a multifunctional cytochrome P450 in the biosynthesis of the cyanogenic glucoside dhurrin. Plant Molecular Biology 36:393−405 doi: 10.1023/A:1005915507497

    CrossRef   Google Scholar

    [40] Pičmanová M, Neilson EH, Motawia MS, Olsen CE, Agerbirk N, et al. 2015. A recycling pathway for cyanogenic glycosides evidenced by the comparative metabolic profiling in three cyanogenic plant species. Biochemical Journal 469:375−89 doi: 10.1042/BJ20150390

    CrossRef   Google Scholar

    [41] Bjarnholt N, Neilson EH, Crocoll C, Jørgensen K, Motawia MS, et al. 2018. Glutathione transferases catalyze recycling of auto-toxic cyanogenic glucosides in sorghum. The Plant Journal 94:1109−25 doi: 10.1111/tpj.13923

    CrossRef   Google Scholar

    [42] Gleadow RM, Møller BL. 2014. Cyanogenic glycosides: synthesis, physiology, and phenotypic plasticity. Annual Review of Plant Biology 65:155−85 doi: 10.1146/annurev-arplant-050213-040027

    CrossRef   Google Scholar

    [43] Burke JJ, Chen J, Burow G, Mechref Y, Rosenow D, et al. 2013. Leaf dhurrin content is a quantitative measure of the level of pre-and postflowering drought tolerance in sorghum. Crop Science 53:1056−65 doi: 10.2135/cropsci2012.09.0520

    CrossRef   Google Scholar

    [44] Emendack Y, Burke J, Laza H, Sanchez J, Hayes C. 2018. Abiotic stress effects on sorghum leaf dhurrin and soluble sugar contents throughout plant development. Crop Science 58:1706−16 doi: 10.2135/cropsci2018.01.0059

    CrossRef   Google Scholar

    [45] Hoagland DR, Arnon DI. 1938. The water culture method for growing plants with soil. Rep. of the Calif. California Agricultural Experiment Station Circular 1950:39

    Google Scholar

    [46] Hiscox JD, Israelstam GF. 1979. A method for the extraction of chlorophyll from leaf tissue without maceration. Canadian Journal of Botany 57:1332−34 doi: 10.1139/b79-163

    CrossRef   Google Scholar

    [47] Arnon DI. 1949. Copper enzymes in isolated chloroplasts. Polyphenoloxidase in Beta vulgaris. Plant Physiology 24:1 doi: 10.1104/pp.24.1.1

    CrossRef   Google Scholar

    [48] Blum A, Ebercon A. 1981. Cell membrane stability as a measure of drought and heat tolerance in wheat. Crop Science 21:43−47 doi: 10.2135/cropsci1981.0011183X002100010013x

    CrossRef   Google Scholar

    [49] Xu Y, Huang B. 2018. Comparative transcriptomic analysis reveals common molecular factors responsive to heat and drought stress in Agrostis stolonifera. Scientific Reports 8:15181 doi: 10.1038/s41598-018-33597-3

    CrossRef   Google Scholar

    [50] Simão FA, Waterhouse RM, Ioannidis P, Kriventseva EV, Zdobnov EM. 2015. BUSCO: assessing genome assembly and annotation completeness with single-copy orthologs. Bioinformatics 31:3210−2 doi: 10.1093/bioinformatics/btv351

    CrossRef   Google Scholar

    [51] Robinson MD, McCarthy DJ, Smyth GK. 2010. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics 26:139−40 doi: 10.1093/bioinformatics/btp616

    CrossRef   Google Scholar

    [52] Young MD, Wakefield MJ, Smyth GK, Oshlack A. 2010. Gene ontology analysis for RNA-seq: accounting for selection bias. Genome Biology 11:R14 doi: 10.1186/gb-2010-11-2-r14

    CrossRef   Google Scholar

    [53] Yu G, Wang L, Han Y, He Q. 2012. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS: A Journal of Integrative Biology 16:284−87 doi: 10.1089/omi.2011.0118

    CrossRef   Google Scholar

    [54] Langfelder P, Horvath S. 2008. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics 9:559 doi: 10.1186/1471-2105-9-559

    CrossRef   Google Scholar

    [55] Zhang B, Horvath S. 2005. A general framework for weighted gene co-expression network analysis. Statistical Applications in Genetics and Molecular Biology 4:17 doi: 10.2202/1544-6115.1128

    CrossRef   Google Scholar

  • Cite this article

    Xu Y, Rossi S, Huang B. 2021. Comparative transcriptomics and gene network analysis revealed secondary metabolism as preeminent metabolic pathways for heat tolerance in hard fescue. Grass Research 1: 12 doi: 10.48130/GR-2021-0012
    Xu Y, Rossi S, Huang B. 2021. Comparative transcriptomics and gene network analysis revealed secondary metabolism as preeminent metabolic pathways for heat tolerance in hard fescue. Grass Research 1: 12 doi: 10.48130/GR-2021-0012
    shu

Figures(8)

Download PDF

Article Metrics

Article views(5174) PDF downloads(815)

Other Articles By Authors

ARTICLE   Open Access    

Comparative transcriptomics and gene network analysis revealed secondary metabolism as preeminent metabolic pathways for heat tolerance in hard fescue

  • Department of Plant Biology, Rutgers, The State University of New Jersey, New Brunswick, NJ 08901, USA

  • Received: 08 October 2021
  • Accepted: 01 December 2021
  • Published online: 30 December 2021
Grass Research  1,  Article number: 12  (2021)  |  Cite this article

Abstract: The molecular mechanisms underlying genetic variations in heat tolerance, one of the important turfgrass traits, for fine fescue are not well-understood. In the present study, our objective was to identify molecular constituents and metabolic interactions involved in heat tolerance in two genotypes of hard fescue (Festuca trachyphylla) contrasting in heat tolerance by comparative transcriptomics and gene comparison network analysis. Two cultivars of hard fescue, 'Reliant IV' (heat-tolerant), and 'Predator' (heat-sensitive), were subjected to heat stress temperature at 35/30 °C (day/night) or maintained at optimal temperature at 22/18 °C (day/night) (non-stress control) for 21 d. At 14 and 21 d of heat stress, 'Reliant IV' maintained significantly higher photochemical efficiency (Fv/Fm), chlorophyll (Chl) content, and lower cell membrane electrolyte leakage (EL) compared to 'Predator', suggesting its superiority in heat tolerance. Comparative transcriptomic profiles, gene functional enrichment analysis, and weighted gene comparison network analysis revealed central hub genes (BBE22 and ALPLD) and their connecting genes involved in secondary metabolism for biosynthesis of oxylipins (LOX1 and LOX3), phenolic compounds (PAL2), and dhurrin (C79A1 and C71E1). These genes were up-regulated in heat-tolerant 'Reliant IV' under heat stress but not in heat-sensitive 'Predator', while a majority of heat-regulated genes involved in primary metabolism responded similarly to heat stress in both cultivars. Those unique genes in the secondary metabolic pathways enriched in only the heat-tolerant cultivar could be critical for mediating the protection of hard fescue against heat stress and are potentially useful as candidate genes or molecular markers for augmenting heat tolerance in other temperate species of grass.

    HTML

INTRODUCTION

  • The growth and productivity of temperate plants is commonly limited by heat stress. Over the next hundred years, the global temperature is projected to rise by 6.9 °C, which could diminish crop yield by 15%−35%[1,2]. Extensive efforts have been made to study the physiological and molecular factors regulating tolerance to heat stress in a multitude of plants (for reviews, see[3−8]). In the past decade, automated equipment designed for rapid testing has become prevalent, generating an appeal for conducting assays seeking to analyze how '–omic' technologies (associated with genes, mRNA, proteins, and metabolites) factor into the heat tolerance of plants.

    Transcriptomic analysis is a powerful approach for the identification of genes associated with plant responses or resistance to abiotic stresses by RNA sequencing (RNA-seq), a technique previously utilized to study responses to heat stress in a variety of temperate grasses, including Kentucky bluegrass (Poa pratensis)[9], perennial ryegrass (Lolium perenne)[10], creeping bentgrass (Agrostis stolonifera)[11], and tall fescue (Festuca arundinacea Schreb.)[12]. Most of the previous studies on transcriptomic analysis in species of cool-season grasses focused on examining transcript changes in a single genotype responding to heat stress which reflected the heat responses of transcriptomes, and a majority of the heat-responsive genes in these studies was associated with primary metabolism of photosynthesis, respiration, proteins, and hormones, as well as transcriptional factors associated with response to heat stress, namely heat shock factors[9−10, 12]. Secondary metabolism for biosynthesis of compounds, such as phenolics, is well known to play roles in plant defense against biotic stresses, but limited information is available on key genes or networks involved in the metabolic pathways of secondary metabolites for plant tolerance to heat stress and are not well documented[3, 5]. Comparative analysis of transcriptomic changes under prolonged periods of heat stress for genotypes differing in tolerance to heat stress, particularly in stress-tolerant grass species, will allow for the classification of unique genes that correlate to genetic variations in heat tolerance, which could be useful for the utilization of genes regulating plant heat tolerance through genomic manipulation or molecular breeding.

    Fine fescues, including hard fescue, are low-maintenance cool-season grass species in terms of fertility, irrigation, and mowing, and also exhibit superb heat tolerance with a broad variety of genetic variability among genotypes. Previous studies[13,14] that compared heat tolerance for 26 genotypes, including seven hard fescues, two slender creeping red fescues (Festuca rubra subsp. litoralis), two sheep fescues (Festuca ovina subsp. hirtula), seven strong creeping red fescues (Festuca rubra subsp. rubra), and eight chewings fescues (Festuca rubra subsp. commutata), demonstrated significant genetic variability in heat tolerance, with 'Reliant IV' hard fescue ranked highly as a heat-tolerant genotype and 'Predator' being heat-sensitive relative to 'Reliant IV', based on the evaluation of turf quality (TQ), photochemical efficiency (Fv/Fm) chlorophyll (Chl) content, relative water content (RWC), and electrolyte leakage (EL). Quantitative PCR analysis showed that there was a direct correlation between the abundance of transcripts for genes associated with energy production (ATP synthase), cellular respiration (Sucrose synthase), light-harvesting (Photosystem II CP47 reaction center protein, RuBisCO activase), oxidative defense (Catalase), cell structure and division (Actin), and stress defense (Heat shock protein 90) and heat tolerance, as manifested by the aforementioned physiological traits[14]. Therefore, we hypothesized that a more-detailed comparative analysis on transcriptomes of such two hard fescue cultivars can reveal a more comprehensive gene network and potential central hub genes that are responsible for their contrasting heat tolerances. For further understanding of the molecular factors and gene networks linked to heat tolerance in hard fescue, we used comparative transcriptomics by RNA-seq and conducted a weighted gene comparison network analysis for the heat-tolerant 'Reliant IV' and heat-sensitive 'Predator' cultivars to establish central hub and unique genes and underlying metabolic pathways that may relate to genetic differences in heat tolerance for hard fescue. Such information will be substantial in the efforts towards genetic augmentation of heat tolerance in cool-season grasses.

    RESULTS

      Differential physiological responses of 'Predator' and 'Reliant IV' to heat stress

    • Under non-stress temperature (0 d of heat stress) and at 7 d of heat stress, 'Predator' and 'Reliant IV' did not differ significantly in canopy density and color (Fig. 1a & b). 'Reliant IV' had more green leaves and a lower proportion of yellow leaves at 14 d of heat stress (Fig. 1c).

      Figure 1.  Plants of 'Predator' and 'Reliant IV' grown under (a) normal temperature (0 d of heat stress), (b) 7 d of heat stress, (c) 14 d of heat stress.

      Leaf Chl content in 'Predator' and 'Reliant IV' decreased with heat stress at 14 and 21 d, and Chl content was significantly greater in 'Reliant IV' as compared with 'Predator' at 14 d of heat stress (Fig. 2a). Leaf Fv/Fm in 'Predator' and 'Reliant IV' declined during heat stress, but 'Reliant IV' maintained significantly higher Fv/Fm than 'Predator' from 7−21 d of heat stress (Fig. 2b). Heat stress caused increases in leaf EL in both 'Predator' and 'Reliant IV', but 'Reliant IV' had significantly lower EL than 'Predator' from 7−21 d of heat stress (Fig. 2c).

      Figure 2.  (a) Leaf chlorophyll (Chl) content, (b) photochemical efficiency (Fv/Fm), and (c) electrolyte leakage (EL) of 'Predator' and 'Reliant IV' in response to heat stress. The error bars represent standard error (SE) for three biological replicates (n = 3).

      Differential transcriptomic responses of 'Predator' and 'Reliant IV' to heat stress

    • The RNA sequencing analysis yielded a total of more than 574 and 513 million reads in 'Predator' and 'Reliant IV', respectively, which was about 90% of the Q30 rate. Both cultivars had over a 97% alignment rate with a similar total contig number, N50 value, average contig length, total assembled bases, GC content (Supplemental Table S1), and gene annotations (Supplemental Table S2), indicating that the transcriptome assembly could mostly represent their transcript profile.

      In 'Predator, a total of 2,184 and 4,324 genes were up-regulated at 7 d while 2,093 and 2,352 genes were down-regulated at 7 d and 14 d of heat stress, respectively (Fig. 3a & b; Supplemental Table S3). Among the differentially expressed genes (DEGs) in 'Predator', 1,792 were commonly up-regulated while 1,585 were commonly down-regulated in response to heat stress at both 7 and 14 d.

      Figure 3.  Venn diagram illustrating the number of differentially expressed genes (DEGs) that were (a) up-regulated or, (b) down-regulation at 7 d or 14 d of heat stress in comparison to that at 0 d of heat stress (Day 7 vs. Day 0 or Day 14 vs. Day 0) for 'Predator'. The numbers in each circle and overlapping area represent unique and common DEGs between days of heat stress, respectively.

      At 7 and 14 d of heat stress, 2,351 and 1,998 genes were up-regulated and 2,227 and 853 were down-regulated, respectively, in 'Reliant IV' (Fig. 4a & b; Supplemental Table S4). Among the DEGs regulated in response to heat stress in 'Reliant IV', 779 were commonly up-regulated and 665 were commonly down-regulated at both 7 and 14 d of heat stress. While a similar number of DEGs were found in 'Reliant IV' and 'Predator' at 7 d of heat stress, fewer up- or down-regulated DEGs were found in 'Reliant IV' than 'Predator' in response to heat stress at 14 d.

      Figure 4.  Venn diagram illustrating the number of differentially expressed genes (DEGs) that were (a) up-regulated or, (b) down-regulation at 7 d or 14 d of heat stress in comparison to that at 0 d of heat stress (Day 7 vs. Day 0 or Day 14 vs. Day 0) for 'Reliant IV'. The numbers in each circle and overlapping area represent unique and common DEGs between days of heat stress, respectively.

      In order to identify major DEGs and their associated biological pathways in relation to genotypic variations in heat tolerance, all up- and down-regulated DEGs at 7 and 14 d of heat stress in 'Predator' and 'Reliant IV' were analyzed using a GO term enrichment program. Among the DEGs up-regulated at 7 d of heat stress, the GO biological process term 'response to oxygen-containing compounds' and GO molecular function term 'oxidoreductase activity' were uniquely enriched in 'Reliant IV', while most enriched GO terms were commonly found in both cultivars, including 'unfolded protein folding', 'heat shock protein folding', 'response to temperature stimulus', 'response to heat', and 'extracellular region' (Fig. 5a; Supplemental Table S5). At 14 d of heat stress, 'Reliant IV' had several uniquely up-regulated and enriched DEGs in 12 GO cellular components, biological processes and molecular terms that were not present in 'Predator', including 'response to oxygen-containing compounds', 'oxidation-reduction process', 'oxylipin biosynthetic process', 'oxylipin metabolic process', 'response to chemical', 'extracellular region', 'intracellular', 'vacuolar lumen', 'ammonia-lyase activity', 'nucleic acid binding transcription factor activity', 'transcriptional factor activity sequence-specific DNA binding', and 'oxidoreductase activity' (Fig. 5b; Supplemental Table S5).

      Figure 5.  Gene ontology (GO) enrichment analysis of DEGs that were up-regulated at (a) 7 d of heat stress and (b) 14 d of heat stress in 'Predator' and 'Reliant IV'.

      The down-regulated DEGs at 7 d of heat stress mostly overlapped between two cultivars, including those categorized in photosynthesis, energy, photosystems, and chlorophyll binding, while DEGs in the biological process 'generation of precursor metabolites and energy', cellular component 'chloroplast part' and 'plastid part', as well as molecular function 'oxidoreductase activity' and 'carbon-carbon lyase activity' categories were down-regulated only in 'Predator' (Fig. 6a; Supplemental Table S5). At 14 d of heat stress, DEGs in four molecular function and biological process categories, including 'photosynthesis', 'ribulose-bisphosphate carboxylase activity', and 'lyase activity', were uniquely down-regulated in 'Predator', while most of the other processes were shared by both cultivars (Fig. 6b, Supplemental Table S5).

      Figure 6.  Gene ontology (GO) enrichment analysis of DEGs that were down-regulated at (a) 7 d of heat stress and (b) 14 d of heat stress in 'Predator' and 'Reliant IV'.

      The KEGG enrichment analysis demonstrated that most of the enriched KEGG terms were identical between 'Predator' and 'Reliant IV' (Fig. 7a & b), whereas the main difference in gene functional enrichment between two cultivars was the up-regulated DEGs at 7 d of heat stress, with unique enrichment in 'Reliant IV' for genes involved in primary and secondary metabolism.

      Figure 7.  KEGG enrichment analysis of DEGs for (a) 'Predator' and (b) 'Reliant IV' at 7 and 14 d of heat stress in comparison to 0 d of heat stress.

      To further explore gene networks and central hub genes related to heat stress tolerance in hard fescue, weighted gene co-expression network analysis (WGCNA) was performed on all heat-responsive DEGs in 'Predator' and 'Reliant IV' exposed to 7 or 14 d of heat stress. The WGCNA grouped DEGs having similar expression arrays under heat stress for 'Predator' and 'Reliant IV' into 15 and 13 module eigengenes, respectively (Supplemental Fig. S1). The gene networks and central hub genes linked to the DEGs up-regulated in only 'Reliant IV' or to DEGs regulated more in 'Reliant IV' with respect to 'Predator' at 7 and 14 d of heat stress, which were enriched in the GO-term functional categories (Fig. 5) for 'response to oxygen-containing compounds', 'oxidation-reduction process', 'oxidoreductase activity', 'oxylipin biosynthetic process', 'oxylipin metabolic process', 'extracellular region', and 'ammonia-lyase activity', were found mainly in two WGCNA modules, MElghtcyan and MEdarkmagenta.

      The central hub genes and connecting genes in the MElightcyan and MEdarkmagenta modules had a log2 fold change in gene expression levels ranging from 1.59 to 7.2 in 'Reliant IV' subjected to heat stress for 7 and 14 d in comparison to that at 0 d of heat stress. The central hub genes in MElightcyan and MEdarkmagenta were BBE22 (Berberine bridge enzyme-like 22 protein) and ALPL (transposase-derived ALP1-like protein), respectively (Fig. 8a & b; Supplemental Table S6 & S7). The unique DEGs in heat-stressed 'Reliant IV' in MElightcyan included LOX1 (lipoxygenase 1) and LOX3 (lipoxygenase 3) in the pathways of oxylipin metabolism, and PAL2 (phenylalanine ammonia-lyase 2) in the processes involving ammonia-lyase activity. The unique DEGs found in MEdarkmagenta included two genes associated with the oxidation-reduction process or oxidoreductase activity for dhurrin biosynthesis, C71E1 (4-hydroxyphenylacetaldehyde oxime monooxygenase) and C79A1 (tyrosine N-monooxygenase), respectively.

      Figure 8.  The heat-upregulated central hub genes and the connecting genes in module (a) MElightcyan and (b) MEdarkmagenta in 'Reliant IV' exposed to 7 and 14 d of heat stress. The genes in the yellow circles represent central hub genes.

    DISCUSSION

    • The correlative analysis of transcriptional alterations as a result of heat stress for two hard fescue cultivars differing in heat tolerance allowed for the identification of several unique heat-responsive DEGs in heat-tolerant 'Reliant IV', while a majority of the heat-responsive DEGs were common in both cultivars, including those with putative functions in primary metabolism, such as photosynthesis and respiration, hormone metabolism, and transcription factors relative to stress that were previously reported in transcriptomic profiling of heat-responsive genes in other plant species[9−12,14]. The WGCNA and GO-term enrichment analysis revealed three unique central hub genes (BBE22 and ALPL) involved in secondary metabolism for stress protection, only up-regulated in 'Reliant IV' under heat stress in this study; however, the manners in which those genes are related to heat tolerance have not previously been examined. Therefore, the discussion below focuses only on those unique heat-responsive DEGs with functions categorized in secondary metabolism, as found in the GO enrichment analysis and the central hubs (BBE22 and ALPL) which may have critical functions in conferring heat tolerance to the heat-tolerant cultivar of hard fescue in order to further understand novel mechanisms of heat tolerance in temperate grasses.

      A central hub gene, BBE22, and three genes, LOX1, LOX3, and PAL, linked with BBE2, were found to be uniquely enriched and up-regulated in the heat-tolerant cultivar, 'Reliant IV', but not in the heat-sensitive 'Predator', under heat stress. The four genes have been associated with secondary metabolism for stress defenses, particularly biotic stress, while their functions for heat stress tolerance are not well-known. LOXs are a family of genes encoding enzymes responsible for catalyzing the hydroperoxidation of polyunsaturated fatty acids, acting on linolenic and linoleic acid in secondary metabolic pathways to synthesize oxylipins, such as jasmonic acid (JA) and green leaf volatiles (GLVs)[15,16]. LOXs play important roles in secondary metabolism, response to naturally-occurring and stress-induced senescence, and defense responses to pathogen and insect attack[17]. LOX1 is inducible by abscisic acid, pathogens, JA[18], and cadmium exposure in Arabidopsis (Arabidopsis thaliana)[19]. Arabidopsis lox1 mutants were shown to have enhanced oxidative damage induced by cadmium[19]. LOX3 encodes an enzyme participating in the synthesis of JA[20] and has previously been exhibited to play a role in JA-mediated salt stress response, as lox3 mutants of Arabidopsis with lower JA content were hypersensitive to salt stress and had slower seed germination rate and lower seedling survival rate, whereas JA treatment rescued the salt-sensitive phenotypes in lox3 mutants[21]. Some studies have also reported that LOX genes (LOX1, LOX3, and LOX5) regulate organ development (i.e. seed germination and lateral rooting) and the levels of LOX activity have been positively associated with cell elongation in rapidly growing tissues of plants[17,22−24]. At the start of the phenylpropanoid pathway, phenylalanine ammonia-lyase (PAL) catalyzes the removal of an amino group from phenylalanine to produce precursors for many secondary metabolites having critical functions in plant defense against biotic and abiotic stresses[25,26]. The PAL2 gene mediates the biosynthesis of lignin and phenolics and was previously reported to be up-regulated in Achillea pachycephala[27] subjected to drought stress and in Arabidopsis[28] exposed to salt stress. The berberine bridge enzyme (BBE) has been associated with plant immunity and nitric oxide responses[29−32]. In a recent study of Arabidopsis response to spider mite, BBE22 was found to be involved in H2O2 homeostasis[33]. lox1 mutants of Arabidopsis exhibited down-regulation of transcription levels of H2O2-scavenging enzymes and significant increases in H2O2 accumulation in leaves and roots exposed to cadmium[19]. The up-regulation of PAL transcripts was associated with increases in the enzymatic activity of peroxidase in barley (Hordeum vulgare L.) exposed to salt stress[34]. The unique up-regulation of the central hub gene of BBE22 linked with LOX1, LOX3, and PAL in the heat-tolerant cultivar of hard fescue suggested that stress defenses involving secondary metabolism could be activated through the central control of BBE22, involving H2O2 scavenging and homeostasis; however, this notation deserves further investigation.

      The ALPL (transposase-derived ALP1-like protein) was another central hub gene upregulated in heat-tolerant 'Reliant IV' exposed to heat stress, which was linked with two DEGs, CYP71E1 (4-hydroxyphenylacetaldehyde oxime monooxygenase) and CYP79A1 (tyrosine N-monooxygenase) in the pathway of dhurrin biosynthesis. Little is known about the biological function of ALPL, except that ALPL1 was identified as a stress-responsive transcription factor related to herbicide resistance[35,36]. The two DEGs linked to ALPL, CYP79A1 and CYP71E1, are localized in the membrane and activate the conversion of tyrosine to dhurrin, which can account for 30% of the dry mass in sorghum (Sorghum bicolor) seedlings[37−39]. Dhurrin can also serve as a stored form of nitrogen that can be remobilized in order to avoid hydrogen cyanide formation and provide a vital source of reduced nitrogen[40,41]. Dhurrin is an insect-defense compound and can be hydrolyzed by β-glucosidase to release toxic hydrogen cyanide gas upon tissue disruption caused by feeding insect larvae[42]. It was previously shown that an accumulation of dhurrin level in leaves is directly correlated to better tolerance to drought in post-flowering sorghum lines[43,44]. The unique upregulation of ALPL, CYP79A1, and CYP71E1 in heat-tolerant cultivars of hard fescue suggests their possible involvement in the heat tolerance of plants; however, the direct link between ALPL, CYP79A1, and CYP71E1 and heat stress tolerance deserves further investigation.

      In summary, unique heat-responsive DEGs found in 'Reliant IV' by functional enrichment and WGCNA analysis include central hub genes (BBE22 and ALPL) and their connecting genes involved in secondary metabolism for biosynthesis of oxylipins (LOX1 and LOX3), phenolic compounds (PAL2), and dhurrin (C79A1 and C71E1), which could be involved in mediating defense against heat stress in plants. The direct causal relationship of those unique heat-responsive genes to heat tolerance warrants further investigation, as they could be useful for the determination of molecular markers associated with superior heat tolerance in hard fescue and other cool-season turfgrass species.

    MATERIALS AND METHODS

      Plant materials and growth conditions

    • Plants of two cultivars of hard fescue, 'Predator' and 'Reliant IV', were transported from a Rutgers University turfgrass research farm (North Brunswick, NJ) to a greenhouse maintained at a daily temperature of 23/20 °C (day/night), with natural and artificial sodium-vapor lighting of 750 μmol m−2 s−1 photosynthetically active radiation (PAR). Experimental plants were transplanted into 15 cm in diameter × 20 cm deep plastic pots filled with an even soil and peat moss mixture. During the 6-week establishment period, plants were irrigated daily, trimmed biweekly to a canopy height of 7 cm, and fertilized biweekly with an even volume of of half-strength Hoagland's nutrient solution[45]. Prior to induction of heat stress treatment, mature plants were transported to controlled-environment growth chambers (Environmental Growth Chambers, Chagrin Falls, OH) set to a 14-h photoperiod at 700 μmol m−2 s−1 PAR, 50% relative humidity, and 22/18 °C (day/night) for plant establishment.

      Temperature treatments and experimental design

    • Plants subjected to heat stress were maintained at a temperature of 38/33 °C (day/night) for 21 d, and non-stress control plants were kept at an optimal temperature of 22/18 °C (day/night). Both temperature regiments were replicated across three growth chambers, and for each cultivar, there were three individual replicates randomly dispersed among the three chambers for each temperature treatment. Plants were irrigated each day and half-strength Hoagland's nutrient solution was applied as a fertilizer weekly. This experiment was arranged as a split-plot randomized design, having temperature conditions defined as the main plots and cultivars defined as the sub-plots.

      Physiological analyses

    • To measure the chlorophyll (Chl) content of leaf tissue, the methods of Hiscox and Israelstam[46] were followed with minor changes made to the procedure. Leaf tissue (100 mg) was excised from plants and submerged in 10 mL dimethyl sulfoxide for 5 d until Chl was completely extracted. After measuring the absorption of the Chl extract at wavelengths of 663 and 645 nm via spectrophotometer (Spectronoic Instruments, Rochester, NY), the leaf tissue was oven-dried for 3 d at 80 °C, and leaf Chl content was determined using the equations provided by Arnon[47].

      A chlorophyll fluorescence meter (FIM 1500, ADC BioScientific Ltd., Herts, UK) was utilized to measure the variable fluorescence (Fv) and maximal fluorescence (Fm) of plants after leaves were dark-adapted for 30 min, and photochemical efficiency (Fv/Fm) was expressed as the ratio of Fv to Fm.

      Leaf electrolyte leakage (EL) was quantified as a measure of membrane stability using the procedure detailed by Blum and Ebercon[48] with certain modifications. Leaf tissue (100 mg) was excised from each plant, submerged in 30 mL of pure, deionized water, and kept on an orbital platform shaker for 16 h. A conductance meter (Model 32, YSI Inc., Yellow Springs, OH) equipped with a probe was used to measure initial conductance (Ci) of each solution. After killing leaf tissue in an autoclave set to a cycle of 20 min at 120 °C, tissue solutions were shaken for 16 h and maximal conductance (Cmax) was measured. EL was calculated by dividing Ci by Cmax and multiplying by 100 to express values as a percentage.

      RNA extraction, library preparation, and RNA sequencing

    • Leaf tissue (0.2 g) was excised from plants at 0, 7, and 14 d of heat stress and total RNA was extracted and pooled into a library, and prepared for RNA sequencing using the methods previously described by Xu and Huang with no modifications[49]. RNA was extracted, and the Low Sample (LS) procedure from the Illumina TruSeq RNA Library Prep Kit v2 (Illumina, San Diego, CA) was utilized to construct pooled libraries before sequencing RNA with a MiSeq Reagent Kit v3 cartridge (Illumina, San Diego, CA).

      Read alignment, counting, and identification of differentially expressed genes (DEGs)

    • Raw sequence reads from MiSeq sequencing were analyzed using the direct methods and program settings provided by Xu and Huang[49], by assembling and trimming reads and analyzing predicted transcripts using the 956 plantae Benchmarking Universal Single-Copy Orthologs (BUSCO) profiles[50]. Transcripts were analyzed for differential expression according to the methods specified in Xu and Huang[49], and edgeR was used for estimating genewise dispersion[51]. To obtain differentially expressed genes (DEGs), the ratios of transcript abundances for 7 d of heat stress to control conditions, 14 d of heat stress to control conditions, and 14 d to 7 d of heat stress conditions for each cultivar were filtered by applying a log2 fold-change threshold of > 2 or < −2 and a false discovery rate (FDR) of < 0.001, and coding regions of assembled transcripts were defined according to the methods specified previously in Xu and Huang[49].

      Gene functional enrichment analysis

    • To determine the gene Ontology (GO) enrichment of DEGs, the GOseq package in R was utilized in order to amend gene length bias with transcriptome data[52]. In order for GO terms to be labeled as 'signifcantly enriched in DEGs', a false discovery rate (FDR) of < 0.05 was required.

      The Kyoto Encyclopedia of Genes and Genomes (KEGG) is a database which includes a hierarchical structure of different levels of biological processes, molecular functions, and cellular components collected from all living organisms. For KEGG enrichment analysis, we used the ClusterProfile R package[53]. Since our species are not available for direct KEGG enrichment, the enricher function in the package was called to use customized KEGG terms in both DEGs and annotations and to calculate the over-represented KEGG terms using default parameters. For visualization of the final results, we used the dotplot function in the package and modified the figure with the ggplot2 function to compare cultivars side-by-side.

      Weighted gene co-expression network analysis (WGCNA)

    • To construct a co-expression network of genes, the WGCNA package in R was used[54]. First, DEGs expressed in at least one pair-wise comparison of heat stress with control condition were included for co-expression network construction. Gene expression levels from all time points were initially clustered into multiple groups to calculate the sample height. A soft threshold of 35 was used to generate adjacency matrix as the scale-free topology criterion, and DEGs were further clustered in the above hierarchical structure using topological overlap method[55]. The module in the gene dendrogram was detected by the dynamic tree cut method, using the parameter of mergeCutHeight = 0.1 and minModuleSize = 30. To further create weighted gene co-expression network, gene connectivity was calculated from the edge weight in the topology overlap matrix, which represents the strength of interactions between two genes. Finally, the weights of all edges of a node (gene) in the network were added together to be the level of connectivity, and the nodes with the highest connectivity in each module, calculated using chooseTopHubInEachModule function in the package, were considered as central hub genes for each module.

      Statistical analysis and data accessions

    • -way analysis of variance was implemented to separate the differences between temperature treatments, cultivars, and their interactive effects on different parameters using SAS (SAS Int. Cary, NC). Differences among treatments and cultivar variations were defined using Student's t-test at a probability level of p = 0.05.

      Transcriptome shotgun assembly for 'Predator' and 'Reliant IV' were installed onto the GenBank Transcriptome Shotgun Assembly (TSA) database, using the accession number PRJNA632630. The reads from all replicates of 'Predator' and 'Reliant IV' are stored in GenBank accession SAMN14916095 - SAMN14916112.

      Acknowledgments
      • The authors wish to thank the National Institute of Food and Agriculture, USDA, Specialty Crop Research Initiative for funding (award number 2017-51181-27222).
      Conflict of interest
      • The authors declare that they have no conflict of interest.
      Supplementary information
      • Supplemental Table S1 RNA-seq library and assembly statistics.
      • Supplemental Table S2 Summary of RNA-seq annotations found in seven databases.
      • Supplemental Table S3 Lists of up- and down-regulated DEGs found in ‘Predator’ at 7 and 14 d of heat stress.
      • Supplemental Table S4 Lists of up- and down-regulated DEGs found in ‘Reliant IV’ at 7 and 14 d of heat stress.
      • Supplemental Table S5 GO term enrichment of DEGs in both ‘Predator’ and ‘Reliant IV’ at 7 and 14 d of heat stress.
      • Supplemental Table S6 WGCNA analysis of DEGs showed two modules (MElightcyan and MEdarkmagenta) and pair-wise gene connection weight in each of them.
      • Supplemental Table S7 Gene expression pattern of the central hub genes and connecting genes in MElightcyan and MEdarkmagenta modules, showing as log of fold change(logFC) ± standard error(SE, n = 3). N.S. means ‘not significant’.
      • Supplemental Fig. S1 Cluster dendrograms of WGCNA analysis for DEGs in ‘Predator’ (a)  and ‘Reliant IV’ (b)  exposed to 7 d and 14 d of heat stress. The resulting original and merged module eigengenes after dynamic tree cut method were also shown at the bottom.
      Rights and permissions
      • Copyright: © 2021 by the author(s). Exclusive Licensee Maximum Academic Press, Fayetteville, GA. This article is an open access article distributed under Creative Commons Attribution License (CC BY 4.0), visit https://creativecommons.org/licenses/by/4.0/.
    Figure (8)  References (55)
  • About this article
    Cite this article
    Xu Y, Rossi S, Huang B. 2021. Comparative transcriptomics and gene network analysis revealed secondary metabolism as preeminent metabolic pathways for heat tolerance in hard fescue. Grass Research 1: 12 doi: 10.48130/GR-2021-0012
    Xu Y, Rossi S, Huang B. 2021. Comparative transcriptomics and gene network analysis revealed secondary metabolism as preeminent metabolic pathways for heat tolerance in hard fescue. Grass Research 1: 12 doi: 10.48130/GR-2021-0012
    shu

Catalog

  • Html

  • Introduction
  • Results

  • Differential physiological responses of 'predator' and 'reliant iv' to heat stress

  • Differential transcriptomic responses of 'predator' and 'reliant iv' to heat stress

  • Discussion
  • Materials and methods

  • Plant materials and growth conditions

  • Temperature treatments and experimental design

  • Physiological analyses

  • Rna extraction, library preparation, and rna sequencing

  • Read alignment, counting, and identification of differentially expressed genes (degs)

  • Gene functional enrichment analysis

  • Weighted gene co-expression network analysis (wgcna)

  • Statistical analysis and data accessions

  • Acknowledgments
  • Conflict of interest
  • Supplementary information
  • Rights and permissions
  • About this article
logo

Contents

Policies

Services

Partnerships

© 2025 Maximum Academic Press, except Open Access articles. MAP is a member of ALPSP, COPE, CrossRef, cOAlition S, OASPA and STM.

Privacy Policy  |  Site Map  |  Terms and Conditions

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return