Search
Submit Manuscript
1.6CiteScore
Advanced Search
MAP  >  Grass Research
2021 Volume 1
Article Contents
Abstract
Introduction
Materials and methods
Plant materials and growth conditions
Identification and phylogenetic analysis of pbwoxs
Characterization of pbwox genes and proteins
Rna extraction and gene expression analysis
Statistical analysis
Results
Identification and phylogenetic analysis of pbwoxs
Physicochemical properties and analysis of conserved motifs of pbwoxs
Prediction of cis-acting elements in the pbwox promoters
Diversified expression patterns of pbwoxs among tissues
Expression patterns of pbwoxs during se of p. bournei
Expression profiles of pbwoxs in response to hormone treatment
Discussion
Characteristics and expression of the wox family members in p. bournei
Expression patterns reveal that pbwox participates in se
Response of pbwox genes to various hormones during se
Conclusions
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

    Phoebe bournei, a rare and endangered protected tree species that is unique to China, which produces excellent material and fragrance, can be used for the production of furniture and as an ornamental tree[1] . However, few natural resources of this species are available, and this species undergoes a long juvenile phase. Currently, seed propagation is the main reproduction technique, and unstable yields driven by fruiting characteristics has a large impact on seedling production[2].

    Somatic embryogenesis (SE) is one of the most important techniques for tree breeding programs, but the mechanism underlying SE is poorly understood[3]. In angiosperms, a mature somatic embryo is induced from embryonic calli and subsequently develops into spherical, heart-shaped, torpedo, and cotyledon-producing embryos[4]. Moreover, regulation of the different stages of SE requires specific cell fate changes, and many transcription factors (TFs) are involved in this process. For example, WUSCHEL (WUS), WUSCHEL-related homeobox (WOX), BABY BOOM (BBM), AGAMOUS-like (AGL), LEAFY COTYLEDON (LEC), Receptor-Like Kinase (SERK), and Vmyb Avian Myeloblastosis Viral Oncogene Homolog (MYB) genes function as indispensable regulators transforming nonembryogenic calli cells into embryogenic calli cells or driving changes between the different developmental stages of SE[5−10]. As such, SE requires precise transcriptional regulation. SE regeneration techniques have been determined for P. bournei[11], and high-quality genomic data of this species have been released by our group[12]. However, the transcriptional regulatory mechanism behind the transitions among the different stages of SE in P. bournei remains elusive.

    An increasing number of studies have shown that WOXs are extensively involved in plant organ regeneration, growth and development, stress responses, and other transcriptional regulatory processes, especially those that occur during SE[13−17]. In Arabidopsis, AtWUS was shown to be expressed at the proembryonic (16-cell) stage and is involved in subsequent maturation during SE[18]. AtWOX2, AtWOX8, and AtWOX9 participate in polarity establishment during early embryonic development, and AtWOX2 is expressed in apical cells, while AtWOX8 and AtWOX9 are expressed specifically in basal cells, which are indispensable for the correct establishment of the apical–basal axis[19,20]. Moreover, PaWOX2 was also shown to be highly expressed in embryogenic cells in Picea abies[21], and the overexpression of PpWOX2 was shown to affect related traits of somatic embryos in Pinus pinaster[17]. In Vitis vinifera, VvWOX2 and VvWOX9 are expressed at high levels during SE and can be used as marker genes for SE[22]. Furthermore, MtWOX9-1 was shown to increase the embryogenic capacity of recalcitrant plant species, e.g., Medicago truncatula[23]. These studies have shown that WOXs are crucial during the process of embryonic development or somatic embryo regeneration. Moreover, in woody plant species, global transcriptomic data and expression analysis have resulted in the identification of WOXs expressed during SE in Dimocarpus longan, hybrid sweetgum, and Elaeis guineensis, suggesting that WOXs are functionally conserved in woody plants species[7,9,24]. Based on this, understanding the dynamic relationship between WOXs and SE in P. bournei is helpful for optimizing the somatic embryo regeneration system and creating a large number of clones rapidly.

    Previous studies have shown that overexpression or ectopic expression of embryogenesis-related TFs can induce the SE process. Another way is to apply exogenous phytohormones[25]. Adding exogenous phytohormones to media can affect the morphology and quality of SE in many species[26]. The interactions between phytohormones and TFs has been under increasing scrutiny. Several studies have shown that some TFs, such as LEC2, BBM, and WUS, are regulated by auxin synthesis, transport, and responses during SE[10,27,28]. Correct establishment of the auxin gradient and PIN1-mediated auxin transport were shown to affect the expression level of WUS, which in turn affected the status of the embryonic calli[27]. LEC2 and BBM transcriptionally regulate the endogenous auxin (IAA) biosynthesis-related genes YUCs/TAAs and increase the DR5 auxin response, further maintaining somatic embryo growth[10,28,29]. Furthermore, LEC2 was shown to bind directly to the early embryonic marker genes WOX2 and WOX3, triggering SE[5,25]. Abscisic acid (ABA) is another important hormone involved in SE, especially during embryo maturation. Application of exogenous ABA to the media was shown to induce embryo maturation and prevent early germination in Carica papaya, Pseudotsuga menziesii, and Phoenix dactylifera[30−32]. Methyl jasmonate (MeJA) plays a function similar to that of ABA in promoting mature SE. MeJA functions synergistically with ABA, but the effects of MeJA cannot replace the effects of ABA[33]. In Liriodendron hybrids, MeJA was shown to increase both SE and the maturation rate and decreased the deformation rate[34]. However, studies on the relationships between MeJA and WOXs are lacking. Taken together, these results suggested that TFs and hormones jointly regulate plant SE. In P. bournei, how WOXs respond to hormones during SE has not been thoroughly characterized. So we preliminarily explored the expression patterns of WOX under auxin, ABA, and MeJA treatments.

    In the present study, 15 WOX genes were identified across the P. bournei genome, and their gene structures and protein sequences were characterized. Then, the expression patterns of WOXs among six tissues and at different stages of SE were determined. To elucidate how these WOXs respond to hormones, their expression levels in response to auxin, ABA, and MeJA were analyzed. Our results revealed WOX members in P. bournei and several possible associations between WOXs and plant hormones. The results of this study will provide further insight into the function of WOXs involved in regulating SE in woody plant.

    MATERIALS AND METHODS

    Plant materials and growth conditions

    The half-sibling family of P. bournei designated 'WY1' was cultivated in the greenhouse. The epicotyls, stem tips, roots, stems, and leaves of three-month-old seedlings and embryogenic calli induced from immature embryos were frozen in liquid N2 and used for semiquantitative analysis of PbWOX genes. Growth of embryonic calli was induced in cotyledon-stage embryos of the 'WY1' mother tree and subcultured at 24 °C in the dark, as described in our previous study[11]. Somatic embryos of P. bournei at six developmental stages, including three stages of calli, the globular embryo stage, and immature and mature cotyledon-producing embryo stages, were collected under a stereomicroscope (OLYMPUS, Beijing, China) and then frozen in liquid N2 for RNA extraction. With respect to calli growth in liquid media for hormone treatment, 0.1 g of calli was transferred to liquid media supplemented with 100 μM IAA, ABA, and MeJA for 3, 6, 12, 24, and 48 h. Calli in untreated liquid media were used as controls. Sampling was performed at the same time, and three replicates were included.

    Identification and phylogenetic analysis of PbWOXs

    The PbWOXs were identified by two methods. Firstly, using the hidden Markov model, we downloaded the sequence of the conserved homeobox domain of the WOX (PFAM00046) from an online website (http://pfam.xfam.org), and the hmm search module in HMMER (version 3.1) software was used to search the protein sequences of the P. bournei genome[12]. The threshold was set to < E−20. Secondly, we downloaded 15 AtWOXs proteins sequence from the TAIR database (www.arabidopsis.org), then used them as query sequences to perform the BLASTp search (E-value < 1e-5) with P. bournei protein sequences. By combining the two methods, candidate sequences without a homeobox domain were omitted. The 15 obtained PbWOX protein sequences were subjected to MUSCLE alignment of MEGA (version 7.0) together with the sequences of 15 AtWOXs, 13 OsWOXs, 18 PtWOXs, eight AtriWOXs, seven SmWOXs, three PpaWOXs, one OstuaWOX, and one OstluWOX protein downloaded from the online Plant Transcription Factor Database (PlantTFDB) website (http://planttfdb.cbi.pku.edu.cn). The neighbor-joining (NJ) method with 1000 bootstrap repetitions was subsequently used to construct a phylogenetic tree, and the other parameters were set to their default.

    Characterization of PbWOX genes and proteins

    Exons and introns of individual PbWOXs were visualized via the online software Gene Structure Display Server (GSDS) (version 2.0) (http://gsds.gao-lab.org), and Multiple Em for Motif Elicitation (MEME) (version 5.11) (http://meme-suite.org/) was used to predict the motifs of the PbWOX family proteins. The ProtParam (https://web.expasy.org/protparam/) online website was subsequently used to predict the physicochemical properties of PbWOX family members, such as their number of amino acids, molecular weight, and isoelectric point. ClustalX (version 1.81) was used for multiple sequence alignment to confirm the presence of WUS-box domain and the homeobox domain. The genome sequence and gene annotation information file was added to the TBtools GFF3 Sequence Extractor submenu, the upstream bases was set to 2000, the upstream CDS 2.0 kb of all genes in P. bournei were obtained. Then 2.0 kb upstream promoter sequences of 15 PbWOX genes were obtained from the TBtools quick fasta extractor submenu[12,35]. Finally, we uploaded the obtained file to an online site PlantCARE (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/) to analyze cis-acting elements.

    RNA extraction and gene expression analysis

    Total RNA was extracted using a RNAprep Pure Plant Kit (TIANGEN, Beijing, China). Then, the RNA was quantified by a Nanodrop ND-1000 spectrophotometer and checked according to the A260/280 nm and A260/A230 nm values. Subsequently, cDNA was synthesized using a PrimeScriptTM RT Reagent Kit with gDNA Eraser (Perfect Real Time) (Takara, Dalian, China), and each RNA sample was 2000 ng. The obtained cDNA was subsequently diluted five times for quantitative RT‒PCR.

    Specific primers of the 15 PbWOX genes were designed using the Primer 3 online website (http://bioinfo.ut.ee/primer3-0.4.0/), and the sequences of these primers are listed in Supplemental Table S1. The expression levels of the PbWOX genes were detected via quantitative RT–PCR and a CFX 96-well Real-Time PCR System (Bio-Rad, USA). The qPCR mixture volume was 10 μL, which comprised 5 μL of 2× ChamQTM SYBR qPCR Master Mix, 0.4 μL of cDNA, 0.2 μL of forward primer, 0.2 μL of reverse primer, and 4.2 μL of ddH2O. The PCR was carried out as follows: predegeneration at 95 °C for 1 min, 45 cycles of denaturation at 95 °C for 10 s followed by annealing at 57 °C for 10 s, and extension at 72 °C for 20 s. PbEF1α was used as an internal control, and the relative gene expression levels were calculated according to the 2−ΔΔCᴛ method[36].

    Statistical analysis

    All the treatments were performed at least three times. The data were subjected to ANOVA and Duncan's multiple range test at the 5% significance level via SPSS (version 26.0) software.

    RESULTS

    Identification and phylogenetic analysis of PbWOXs

    After performing hidden Markov model (HMM) searches and removing redundant and/or sequences without the homeobox domain, we identified 15 PbWOX members. Phylogenetic analysis of 15 AtWOXs, 18 PtWOXs, 13 OsWOXs, eight AtriWOXs, seven SmWOXs, three PpaWOXs and WOX protein sequences from two green algal species resulted in the assignment of the 15 PbWOX genes to an ancient branch, an intermediate branch and a modern/WUS branch (Table 1). Specifically, the ancient branch consisted of three PbWOXs (PbWOX13a, PbWOX13b, and PbWOX13c); the intermediate branch consisted of four PbWOXs, namely, PbWOX9, PbWOX11/12a, PbWOX11/12b, and PbWOX11/12c, which were classified into two subclasses; and the remaining eight members, namely PbWUS, PbWOX1a, PbWOX1b, PbWOX2a, PbWOX2b, PbWOX3, PbWOX4 and PbWOX5/7, were assigned to the modern/WUS branch (Fig. 1a).

    Figure 1.  Phylogenetic relationships of PbWOX proteins. (a) NJ tree constructed of the amino acid sequence of WOXs from Phoebe bournei (Pb), Arabidopsis thaliana (At), Populus trichocarpa (Pt), Oryza sativa (Os), Amborella trichopoda (Atri), Selaginella moellendorffii (Sm), Physcomitrella patens (Ppa), Ostreococcus tauri (Ostau) and Ostreococcus lucimarinus (Ostlu). (b) Synteny analysis of WOX genes between P. bournei and A. thaliana. Gray lines indicate all synteny blocks in the genome, and the red lines indicate duplicated WOX gene pairs.
    DownLoad: Full-Size Img PowerPoint
    Table 1.  Subclass information of WOXs among P. bournei and other representative species.
    Taxonomic groupSpeciesAncient cladeIntermediate cladeModern/ WUS
    clade    		Total
    DicotsA. thaliana34815
    P. trichocarpa671118
    MonocotsO. sativa16613
    MagnolialesP. bournei34815
    AmborellalesA. trichopoda1258
    PteridophytaS. moellendorffii61−7
    BryophytaP. patens3−−3
    ChlorophytaO. tauri1−−1
    O. lucimarinus1−−1
     | Show Table
    DownLoad: CSV

    However, the number of PbWOX genes was the same as that of Arabidopsis (Fig. 1a). Nonetheless, PbWOXs probably expanded differently than did those of Arabidopsis. For example, three homologs of AtWOX11/12 and AtWOX13, two homologs of AtWOX1 and AtWOX2, one homolog each of AtWUS, AtWOX3, AtWOX4, AtWOX5/7, and AtWOX9, and no homologs of AtWOX6 and AtWOX10 were found in P. bournei (Fig. 1b).

    Physicochemical properties and analysis of conserved motifs of PbWOXs

    A sequence analysis of the PbWOXs showed that PbWOX1b comprised the largest number of amino acid residues (528) and had the largest molecular weight (59.19 kD). Conversely, PbWOX5/7 comprised 169 amino acid residues and had the smallest molecular weight (19.37 kD). All PbWOX genes contain introns, the number of which ranged from two to eight (Fig. 2a). Then, to better understand each member of the PbWOXs, we predicted the physicochemical properties by the use of an online website. The theoretical isoelectric point of PbWOX was found to be between 5.48 (PbWOX11/12c) and 9.93 (PbWOX13c) (Table 2).

    Figure 2.  Information on the PbWOX genes and proteins. (a) Phylogenetic tree and gene structure. (b) Architecture of conserved protein motifs. (c) Multiple sequence alignment.
    DownLoad: Full-Size Img PowerPoint
    Table 2.  Summary of the PbWOX gene family members.
    Gene IDGene nameOrthologous in ArabidopsisTheoretical pIMolecular weightNumber of amino acids
    OF24054-RAPbWUSAtWUS8.5831622.96276
    OF03970-RAPbWOX1aAtWOX19.3737209.6328
    OF11837-RAPbWOX1bAtWOX18.8959188.97528
    OF19048-RAPbWOX2aAtWOX27.0924524.5218
    OF05256-RAPbWOX2bAtWOX26.8324496.61219
    OF16243-RAPbWOX3AtWOX39.0522752.66194
    OF04424-RAPbWOX4AtWOX48.2524797.85220
    OF05362-RAPbWOX5/7AtWOX5, AtWOX79.5119371.68169
    OF24594-RAPbWOX9AtWOX97.1945272.57413
    OF22069-RAPbWOX11/12aAtWOX11, AtWOX125.6830055.89268
    OF11766-RAPbWOX11/12bAtWOX11, AtWOX125.9530330.35277
    OF28194-RAPbWOX11/12cAtWOX11, AtWOX125.4827450.91252
    OF25757-RAPbWOX13aAtWOX135.9132705.61288
    OF14063-RAPbWOX13bAtWOX136.1032380.27286
    OF07768-RAPbWOX13cAtWOX139.9331294.66272
     | Show Table
    DownLoad: CSV

    Motif 1 and motif 2 were detected in all 15 PbWOXs, motif 3 was specific to the members of the intermediate clade, motif 4 (T-L-X-L-F-P-X-X, where X indicates any amino acid) was present in all members of the modern/WUS clade except PbWOX5/7, and motif 5 was specific to PbWOX13a and PbWOX13b (Fig. 2b). There are residues composing homeobox domain motifs that contain three helixes spaced by one loop and one turn (Fig. 2c). Eight members in the modern/WUS clade shared a WUS-box domain (Fig. 2c).

    Prediction of cis-acting elements in the PbWOX promoters

    The cis-acting elements in the promoter region of PbWOXs were divided into four main categories: light-related, hormone-related, stress-related and development-related. (Fig. 3). Specifically, the light response elements constituted the largest proportion, of which the number of G-box elements was the largest. Several other elements involved circadian rhythm were also detected. The hormone-responsive elements included 45 ABA-responsive elements (ABREs), 30 MeJA-responsive elements (CGTCA motif–containing elements), 24 gibberellin (GA)-responsive elements (P-boxes, GARE motif–containing elements, TATC-boxes), 10 salicylic acid-responsive elements (TCA-elements), and nine auxin-responsive elements (TGA-elements, AuxREs, AuxRR-core elements). Abiotic stress response elements were predicted with 38 regulatory anaerobic inductor elements (ARE), 20 drought-responsive elements that could bind MYBs (MBSs), 15 low-temperature–responsive elements (LTRs), eight defense- and stress-responsive elements and five anoxic-specific induction-responsive elements. Moreover, in development-related cis-acting elements, 14 CAT boxes, 12 O2-sites, and six RY elements were predicted, respectively. In the PbWOX promoters, the most common cis-acting elements were G-boxes (light-related), ABREs (ABA-related), CGTCA motif-containing elements (MeJA-related) and AREs (drought-related). This result implied that PbWOX participated in plant growth process and stress response.

    Figure 3.  Predicted cis-acting elements in PbWOX promoters. (a) Frequency of cis-acting elements in the 2.0 kb upstream regions of PbWOXs. The corresponding colored bar chart indicates the occurrence of different cis-acting elements. (b) Number of cis-acting elements in each WOX gene.
    DownLoad: Full-Size Img PowerPoint

    Diversified expression patterns of PbWOXs among tissues

    To further understand the potential roles of PbWOXs during different developmental stages and at different physiological status, semi-qPCR was used to study the expression patterns of 15 PbWOXs in six tissues. The expression levels of PbWOXs varied significantly among the tissues (Fig. 4). Specifically, five genes, namely, PbWOX2a, PbWOX5/7, PbWOX9, PbWOX13a, and PbWOX13b, were expressed in almost all the tissues, while PbWUS, PbWOX1a, PbWOX2b and PbWOX3 were highly expressed in the epicotyls, with low or no expression in the other tissues. In addition, PbWOX11/12a, PbWOX11/12b and PbWOX11/12c were highly expressed specifically in both the roots and embryogenic calli, while expression of PbWOX1b and PbWOX4 was nearly absent in the calli. In total, nine PbWOXs were expressed in embryogenic calli, and thus, these genes may be involved in the SE of P. bournei; PbWOX2a exhibited the highest expression level.

    Figure 4.  Semiquantitative analysis of PbWOXs in different tissues. (a) Tissue samples, 1 - epicotyl, 2 - stem tip, 3 - root, 4 - stem, 5 - leaf, 6 - calli. (b) Semiquantitative PCR electropherogram.
    DownLoad: Full-Size Img PowerPoint

    Expression patterns of PbWOXs during SE of P. bournei

    Previous studies have shown that WOXs play important roles during SE. The expression levels of nine PbWOXs were analyzed in calli at three different developmental stages (Fig. 5a−c) and in embryos at three different developmental stages (Fig. 5d−g). Embryonic calli were induced by immature zygotic embryos (Fig. 5a); then, the embryonic calli developed to the second stage (Fig. 5b) after two or three rounds of propagation, and the calli developed to the third stage (Fig. 5c) after two rounds of propagation. Globular embryos (Fig. 5d), immature cotyledon-producing embryos (Fig. 5e) and mature cotyledon-producing embryos (Fig. 5f) were also selected. The qPCR results showed that the expression levels of PbWOX2a and PbWOX9 increased during embryogenic calli development but decreased as the embryos matured. PbWUS was specifically and highly expressed in the immature cotyledon-producing embryos. The expression level of PbWOX5/7 increased during calli development but decreased after calli differentiation. Three homologous genes, PbWOX11/12 and PbWOX13a, were highly expressed in cotyledon-producing embryos, and their expression peaked upon maturity (Fig. 5g).

    Figure 5.  Expression patterns of PbWOXs during SE of P. bournei. (a) Calli-1. (b) Calli-2. (c) Calli-3. (d) Globular embryo. (e) Immature cotyledon-producing embryo. (f) Mature cotyledon-producing embryo. (g) Analysis of gene expression via qPCR. The data are the means ± SDs of three biological replicates. The values followed by the same letter are not different according to Duncan’s multiple-range test. PbEF1α was used as an endogenous control.
    DownLoad: Full-Size Img PowerPoint

    Expression profiles of PbWOXs in response to hormone treatment

    With respect to the cis-acting elements of PbWOXs, we investigated the expression patterns of PbWOXs in response to auxin, ABA, and MeJA (Fig. 6). Under IAA treatment, PbWUS expression was induced and increased continuously as the treatment duration increased; PbWOX5/7 was strongly induced after 3 h of treatment, after which the expression level decreased. The expression levels of PbWOX2a and PbWOX9 significantly decreased, and the expression levels of PbWOX11/12b, PbWOX11/12c, PbWOX13a, and PbWOX13b also slightly decreased.

    Figure 6.  Relative expression levels of PbWOXs under hormone treatment. The data are the means ± SDs of three biological replicates. The values followed by the same letter are not different according to Duncan's multiple-range test. PbEF1α was used as an endogenous control.
    DownLoad: Full-Size Img PowerPoint

    PbWUS was also induced in response to ABA treatment, while PbWOX2a, PbWOX9, and PbWOX13b were inhibited. The expression levels of PbWOX5/7 and PbWOX13a decreased, reached their lowest level after 12 h of ABA treatment, and then gradually increased. PbWOX11/12b and PbWOX11/12c showed similar expression patterns; their expression increased after 3 h but then decreased. PbWUS expression was induced in response to MeJA treatment, peaked at 12 h, and then gradually decreased. PbWOX5/7, PbWOX9 and PbWOX13b expression was inhibited significantly. PbWOX11/12b and PbWOX11/12c expression increased after 3 h but then decreased.

    DISCUSSION

    Characteristics and expression of the WOX family members in P. bournei

    WOXs are specific to plants and largely involved in key developmental processes, especially those associated with somatic cell regeneration. With the publication of many plant genome sequences, WOX genes have been identified in several plant species. In the present study, we identified 15 PbWOXs, same as the number in Arabidopsis[37], and different orthologous revealed that chromosomal duplication events may occur in P. bournei. Furthermore, the length of introns showed regular characteristics across different clades. For instance, genes in the intermediate clade contained shorter intron sequences than did those of ancient clade, and five genes in the modern/WUS clade had the shortest intron. Taken together, these results suggested that the intron fragments underwent refinement during the evolution of the PbWOX genes. A similar phenomenon was observed in Camellia sinensis[38], which was exemplified by most WOX introns in members of the modern/WUS clade are much shorter than those in the ancient clade. In addition, compared with that in algae, ferns and other more ancestral plant species with one or two members, the WOX family gene in woody plant species has expanded in number and evolved in terms of sequence.

    Tissue-specific expression of a gene implies that the gene plays an indispensable role in certain tissues. We found that, like those in Arabidopsis, the WUS genes in P. bournei were mainly expressed in the epicotyls and shoot apical meristems (SAMs), but this is unlike the patterns of other popular genes, which are expressed in the SAMs, roots, stems, and leaves[39]. These results suggested that PbWUS might play a crucial role in maintaining the differentiation of the SAM. In Arabidopsis, AtWOX4 participates in TDIF-TDR-WOX4 signaling to maintain vascular meristem organization during secondary growth[40], which is similar to what occurs in poplar[41,42]. Here, PbWOX4 was also highly expressed in the stems; thus, this gene may have a function in P. bournei like that of its homologs in Arabidopsis and poplar. In addition, PbWOX4 was also expressed in the roots, leaves, and other plant tissues except embryogenic calli, suggesting that this gene is not involved in plant regeneration or development in vitro. Like AtWOX11, PtoWOX11/12a, and PtoWOX11/12b, three members, namely, PbWOX11/12a, PbWOX11/12b, and PbWOX11/12c, were expressed in the roots and embryogenic calli. Previous studies have shown that WOX11 is involved in adventitious root formation, which has an essential function in root regeneration during de novo plant regeneration[43−46]. Therefore, it was speculated that these three PbWOX11 members might participate in calli propagation and/or root regeneration in P. bournei.

    Expression patterns reveal that PbWOX participates in SE

    SE is one of the important mechanisms of plant asexual reproduction and is subject to complex transcriptional regulation, which in turn enables precise cell fate transitions and the formation of a complete plant. This hierarchical transcriptional regulatory network structure for SE has been revealed in Arabidopsis; in this process, WOX2 and WOX3 are the key TFs that induce SE[5]. According to the tissue expression patterns among tissues, we identified nine WOX genes that were expressed in embryogenic calli—the early stage of SE.

    WUS plays a crucial role in embryogenesis by promoting the fate of cells to transform and develop into embryos, and WUS can also drive the activity of embryonic stem cells. An earlier study showed that WUS is expressed in the four inner apical cells of 16-cell embryos and promotes the formation of the SAM during embryo development, and overexpression of WUS promotes the formation of high-frequency SE. Like in other species, such as Coffea canephora[47], Medicago truncatula[48], and Gossypium hirsutum[15], WUS overexpression resulted in an increased SE induction ratio. In our study, the expression level of WUS significantly increased in the late stage of embryogenic calli and increased significantly again at the immature cotyledon-producing embryo stage. These results indicated that PbWUS promotes the proliferation of embryogenic calli, affects the establishment of cell axial polarity, such as the formation of apical bud meristems during embryonic development in plants, and especially promotes the transition to cotyledon-producing embryos.

    In addition to WUS, WOX2 and WOX9 are the most reported WOX genes involved in plant SE. In the present study, PbWOX2a and PbWOX9 exhibited similar expression patterns, which were exemplified by higher expression levels observed at the embryogenic cell stage and during early somatic embryo formation. In Arabidopsis, it has been proposed that AtWOX2 and AtWOX9 play crucial roles in apical–basal axis formation during embryo development[19]. AtWOX2 is expressed in the apical cell, whereas AtWOX9 is expressed in the basal cell. These genes expressed at specific sites drive the fate of cells in the embryo. In grapevine, both VvWOX2 and VvWOX9 are labeled marker genes of early embryogenic phases[22]. In addition, WOX2 and WOX9 were found to play crucial roles in the early stage of SE in the gymnosperm Picea abies[21,49]. Therefore, PbWOX2a and PbWOX9 might be marker genes for early embryonic development of P. bournei.

    WOX11 has been reported to be an important upstream gene involved in the generation of root system architecture and to promote adventitious root formation during de novo root organogenesis from leaf explants[44,50], but this gene has not been found to be related to root regeneration in plant SE. We noted that PbWOX11/12a, PbWOX11/12b, and PbWOX11/12c were all detected in the embryogenic calli, specifically in immature and mature cotyledon-producing embryos. A similar phenomenon has been observed in grapevine, exemplified by VvWOX11 being highly expressed in torpedo-stage and cotyledon-producing embryos[22]. These findings further support that WOX11/12 might play an important role in the later stage of somatic embryo development and is probably related to root development, but whether WOX11/12 is involved in root primordium formation remains to be confirmed.

    In P. bournei, two WOX13 genes orthologous to AtWOX13 were detected in calli and somatic embryos, but their expression patterns differed. The expression level of PbWOX13a gradually increased with embryo development, while that of PbWOX13b showed no significant change. To our knowledge, WOX13 is expressed ubiquitously and participates in calli formation and organ reconnection in Arabidopsis[51]. However, the molecular regulatory roles of WOX13 during somatic embryo regeneration remain unclear, and expression profiles have been reported in only Vitis vinifera, in which three VvWOX13 genes exhibited low expression levels in somatic embryos, and the expression profile was unaffected by environmental changes[22]. Our data showed that two PbWOX13s also exhibited ubiquitous expression patterns. Nevertheless, PbWOX13b expression seemingly changed nonsignificantly during embryonic calli induction and mature cotyledon-producing embryos, while the expression level of PbWOX13a slightly increased in the later stage of somatic embryo development. Taken together, these results suggested that PbWOX13a might play a regulatory role at the later stage of somatic embryo development.

    Response of PbWOX genes to various hormones during SE

    SE is a highly efficient method for plant regeneration[52]. Overexpression of WUS, WOX2, WOX9, BBM, and SERK is an efficient way to induce SE, and application of plant growth regulators such as auxin, MeJA, ABA, and GA is another useful method[6,47,53−55]. These hormones undergo crosstalk with various TFs and play a primary role in SE[20]. However, information on interactions between phytohormones and WOX genes in P. bournei is lacking. In our study, referring to the information of cis-acting elements in the promoters of PbWOXs, we evaluated that the expression profiles of PbWOXs in embryogenic calli after treatment with IAA, MeJA, and ABA.

    Auxin was first discovered to affect embryonic initiation in carrot and has been widely used to induce SE not only in angiosperms but also in gymnosperms[56]. Moreover, studies have indicated that auxin distribution is positively correlated with the accumulation of WUS, WOX2, and WOX9 transcripts[27,57,58]. In P. bournei, the expression of WUS and WOX5/7 was induced by auxin. However, the expression of WOX2 and WOX9 was inhibited, opposite to what has been reported in Picea abies[49]. In view of this phenomenon, we analyzed the possible causes of species differences or differences in auxin concentration. Whether WOX2/WOX9 and auxin play a synergistic or antagonistic role in somatic embryo initiation remains to be determined.

    ABA involvement in embryo development and maturation has been demonstrated in the SE of several species. In late embryonic development, LEA proteins accumulate in large quantities and act as components of ABA-inducible systems. On the other hand, exogenous ABA in culture media has been shown to promote the maturation and regeneration of somatic embryos[55]. Previous studies have shown that embryo cells cultured in media supplemented with 100 μM ABA produced more embryos in sugi[59]. Six SlWOXs were significantly upregulated after 3 h of 100 μM ABA treatment in tomato[60]. Our data showed that PbWOX11/12a, PbWOX11/12b and PbWOX11/12c were also briefly induced after 3 h of ABA treatment. At the same time, these three genes were highly expressed in the cotyledon-producing embryo stage of SE. It was further speculated that the expression of PbWOX11/12s is likely to be activated by ABA signaling, thereby promoting somatic embryo maturation.

    MeJA is another hormone that increases somatic embryo induction and maturation rate. The effect of MeJA is similar to that of ABA to some extent, but it cannot replace ABA[33]. Previously, 50 μM and 100 μM MeJA were used to treat embryonic calli of longan[61], and exogenous applications of 10-400 μM MeJA produced more mature somatic embryos to different extents[62] Here, we used qPCR to measure the expression changes of PbWOXs after 100 μM MeJA treatment and hoped to determine the relationship between PbWOXs and MeJA. Our data showed that the expression levels of PbWUS and PbWOX11/12a were induced rapidly after MeJA treatment. In addition, the expression levels of PbWOX2a and PbWOX9 were inhibited by MeJA after 3 h. Therefore, applying MeJA to calli for a suitably short time might promote the somatic embryo differentiation process in P. bournei.

    CONCLUSIONS

    The WOX family is unique to plants, and WOX members play important regulatory roles in plant development, such as embryonic patterning. In the present study, we identified 15 PbWOX members in P. bournei, and their expression patterns among different tissues and SE process were determined, and the relationships between PbWOXs and hormones were also analyzed. These results are helpful to further study the regulatory roles of PbWOXs during SE, thus provides the important gene resources for regulating the SE process in P. bournei and other forestry trees.

    ACKNOWLEDGMENTS

    We acknowledge Wenting Xu from Zhejiang A&F University for providing basic experimental materials. We thank for professor Longjun Cheng (Zhejiang A&F University) for his guidance. This work was supported by the Zhejiang Science and Technology Major Program on Agricultural New Variety Breeding (2021C02070-10), the National Natural Science Foundation of China (32171828 and 32101545).

  • 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(5216) PDF downloads(817)

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