ARTICLE   Open Access    

Transcriptomic analysis of zoysiagrass (Zoysia japonica) provides novel insights into the molecular basis of cold acclimation

More Information
  • Zoysia japonica (Z. japonica) is a warm-season perennial turfgrass commonly grown in the southeastern United States for its relatively low input requirements and general tolerance to drought, shade, and salinity. Improving freezing tolerance is critical for Z. japonica, as it could extend the northern boundaries that the species is able to be grown in. To deepen our understanding of the molecular basis of freezing tolerance in Z. japonica, a transcriptomic approach was utilized to identify genes involved in cold acclimation. 'Meyer', a freeze tolerant cultivar, and 'Victoria', a freeze susceptible cultivar were subjected to cold acclimation and non-cold acclimation treatments to determine the number of differentially expressed genes (DEGs). In response to cold acclimation, a total of 4,609 DEGs were up-regulated and 3,605 DEGs were down-regulated in 'Meyer', while 3,758 DEGs were up-regulated and 3,516 DEGs were down-regulated in 'Victoria'. GO and KEGG enrichment analysis showed several diverse pathways and biological processes, including photosynthesis, transmembrane transport, and plant hormone signal transduction. Integrating these information with previous studies on proteomics and QTL mapping, several candidate genes were identified to be associated with cold acclimation and freezing tolerance across different studies, such as LEA, CIPK, POD, HSF, HSP, MPK, PSII and multiple transcription factors. The candidate genes identified in this study showed great value in potentially being a target for future selection efforts for freeze-tolerant Z. japonica.
  • Crops require a variety of nutrients for growth and nitrogen is particularly important. Nitrogen is the primary factor limiting plant growth and yield formation, and it also plays a significant role in improving product quality[14]. Nitrogen accounts for 1%−3% of the dry weight of plants and is a component of many compounds. For example, it is an important part of proteins, a component of nucleic acids, the skeleton of cell membranes, and a constituent of chlorophyll[5,6]. When the plant is deficient in nitrogen, the synthesis process of nitrogen-containing substances such as proteins decrease significantly, cell division and elongation are restricted, and chlorophyll content decreases, and this leads to short and thin plants, small leaves, and pale color[2,7,8]. If nitrogen in the plant is in excess, a large number of carbohydrates will be used for the synthesis of proteins, chlorophyll, and other substances, so that cells are large and thin-walled, and easy to be attacked by pests and diseases. At the same time, the mechanical tissues in the stem are not well developed and are prone to collapse[3,8,9]. Therefore, the development of new crop varieties with both high yields and improved nitrogen use efficiency (NUE) is an urgently needed goal for more sustainable agriculture with minimal nitrogen demand.

    Plants obtain inorganic nitrogen from the soil, mainly in the form of NH4+ and nitrate (NO3)[1013]. Nitrate uptake by plants occurs primarily in aerobic environments[3]. Transmembrane proteins are required for nitrate uptake from the external environment as well as for transport and translocation between cells, tissues, and organs. NITRATE TRANSPORTER PROTEIN 1 (NRT1)/PEPTIDE TRANSPORTER (PTR) family (NPF), NRT2, CHLORIDE CHANNEL (CLC) family, and SLOW ACTIVATING ANION CHANNEL are four protein families involved in nitrate transport[14]. One of the most studied of these is NRT1.1, which has multiple functions[14]. NRT1.1 is a major nitrate sensor, regulating many aspects of nitrate physiology and developmental responses, including regulating the expression levels of nitrate-related genes, modulating root architecture, and alleviating seed dormancy[1518].

    There is mounting evidence that plant growth and development are influenced by interactions across numerous phytohormone signaling pathways, including abscisic acid, gibberellins, growth hormones, and cytokinins[3,19,20]. To increase the effectiveness of plant nitrogen fertilizer application, it may be possible to tweak the signaling mediators or vary the content of certain phytohormones. Since the 1930s, research on the interplay between growth factors and N metabolism has also been conducted[3]. The Indole acetic acid (IAA) level of plant shoots is shown to decrease in early studies due to N shortage, although roots exhibit the reverse tendency[3,21]. In particular, low NO3 levels caused IAA buildup in the roots of Arabidopsis, Glycine max, Triticum aestivum, and Zea mays, indicating that IAA is crucial for conveying the effectiveness of exogenous nitrogen to the root growth response[20,22,23].

    Studies have shown that two families are required to control the expression of auxin-responsive genes: one is the Auxin Response Factor (ARF) and the other is the Aux/IAA repressor family[2426]. As the transcription factor, the ARF protein regulates the expression of auxin response genes by specifically binding to the TGTCNN auxin response element (AuxRE) in promoters of primary or early auxin response genes[27]. Among them, rice OsARF18, as a class of transcriptional repressor, has been involved in the field of nitrogen utilization and yield[23,28]. In rice (Oryza sativa), mutations in rice salt tolerant 1 (rst1), encoding the OsARF18 gene, lead to the loss of its transcriptional repressor activity and up-regulation of OsAS1 expression, which accelerates the assimilation of NH4+ to Asn and thus increases N utilization[28]. In addition, dao mutant plants deterred the conversion of IAA to OxIAA, thus high levels of IAA strongly activates OsARF18, which subsequently represses the expression of OsARF2 and OsSUT1 by directly binding to the AuxRE and SuRE promoter motifs, resulting in the inhibition of carbohydrate partitioning[23]. As a result, rice carrying the dao has low yields.

    Apples (Malus domestica) are used as a commercially important crop because of their high ecological adaptability, high nutritional value, and annual availability of fruit[29]. To ensure high apple yields, growers promote rapid early fruit yield growth by applying nitrogen. However, the over-application of nitrogen fertilizer to apples during cultivation also produces common diseases and the over-application of nitrogen fertilizer is not only a waste of resources but also harmful to the environment[29]. Therefore, it is of great significance to explore efficient nitrogen-regulated genes to understand the uptake and regulation of nitrogen fertilizer in apples, and to provide reasonable guidance for nitrogen application during apple production[30]. In this study, MdARF18 is identified which is a key transcription factor involved in nitrate uptake and transport in apples and MdARF18 reduces NO3 uptake and assimilation. Further analysis suggests that MdRF18 may inhibit the transcriptional level of MdNRT1.1 promoter by directly binding to its TGTCTT target, thus affecting normal plant growth.

    The protein sequence of apple MdARF18 (MD07G1152100) was obtained from The Apple Genome (https://iris.angers.inra.fr/gddh13/). Mutant of arf18 (GABI_699B09) sequence numbers were obtained from the official TAIR website (www.arabidopsis.org). The protein sequences of ARF18 from different species were obtained from the protein sequence of apple MdARF18 on the NCBI website. Using these data, a phylogenetic tree with reasonably close associations was constructed[31].

    Protein structural domain prediction of ARF18 was performed on the SMART website (https://smart.embl.de/). Motif analysis of ARF18 was performed by MEME (https://meme-suite.org/meme/tools/meme). Clustal was used to do multiple sequence comparisons. The first step was accessing the EBI web server through the Clustal Omega channel. The visualization of the results was altered using Jalview, which may be downloaded from www.jalview.org/download.[32]

    The apple 'Orin' callus was transplanted on MS medium containing 1.5 mg·L−1 6-benzylaminopurine (6-BA) and 0.5 mg·L−1 2,4 dichlorophenoxyacetic acid (2,4-D) at 25 °C, in the dark, at 21-d intervals. 'Royal Gala' apple cultivars were cultured in vermiculite and transplanted at 25 °C every 30 d. The Arabidopsis plants used were of the Columbia (Col-0) wild-type variety. Sowing and germinating Arabidopsis seeds on MS nutrient medium, and Arabidopsis seeds were incubated and grown at 25 °C (light/dark cycle of 16 h/8 h)[33].

    The nutrient solution in the base contained 1.0 mM CaCl2, 1.0 mM KH2PO4, 1.0 mM MgSO4, 0.1 mM FeSO4·7H2O 0.1 mM Na2EDTA·2H2O, 50 μM MnSO4·H2O, 50 μM H3BO3, 0.05 μM CuSO4·5H2O, 0.5 μM Na2MoO4·2H2O, 15 μM ZnSO4·7H2O, 2.5 μM KI, and 0.05 μM CoCl·6H2O, and 0.05 μM CoCl·6H2O, and 0.05 μM CoCl· 6H2O. 2H2O, 15 μM ZnSO4·7H2O, 2.5 μM KI and 0.05 μM CoCl·6H2O, and 0.05 μM CoCl·6H2O, supplemented with 0.5 mM, 2 mM, and 10 mM KNO3 as the sole nitrogen source, and added with the relevant concentrations of KCl to maintain the same K concentration[33,34].

    For auxin treatment, 12 uniformly growing apple tissue-cultured seedlings (Malus domestica 'Royal Gala') were selected from each of the control and treatment groups, apple seedlings were incubated in a nutrient solution containing 1.5 mg·L−1 6-BA, 0.2 mg·L−1 naphthalene acetic acid, and IAA (10 μM) for 50 d, and then the physiological data were determined. Apple seedlings were incubated and grown at 25 °C (light/dark cycle of 16 h/8 h).

    For nitrate treatment, Arabidopsis seedlings were transferred into an MS medium (containing different concentrations of KNO3) as soon as they germinated to test root development. Seven-day-old Arabidopsis were transplanted into vermiculite and then treated with a nutrient solution containing different concentrations of KNO3 (0.5, 2, 10 mM) and watered at 10-d intervals. Apple calli were treated with medium containing 1.5 mg·L−1 6-BA, 0.5 mg·L−1 2,4-D, and varying doses of KNO3 (0.5, 2, and 10 mM) for 25 d, and samples were examined for relevant physiological data. Apple calli were subjected to the same treatment for 1 d for GUS staining[35].

    To obtain MdARF18 overexpression materials, the open reading frame (ORF) of MdARF18 was introduced into the pRI-101 vector. To obtain pMdNRT1.1 material, the 2 kb segment located before the transcription start site of MdNRT1.1 was inserted into the pCAMBIA1300 vector. The Agrobacterium tumefaciens LBA4404 strain was cultivated in lysozyme broth (LB) medium supplemented with 50 mg·L−1 kanamycin and 50 mg·L−1 rifampicin. The MdARF18 overexpression vector and the ProMdNRT1.1::GUS vector were introduced into Arabidopsis and apple callus using the flower dip transformation procedure. The third-generation homozygous transgenic Arabidopsis (T3) and transgenic calli were obtained[36]. Information on the relevant primers designed is shown in Supplemental Table S1.

    Plant DNA and RNA were obtained using the Genomic DNA Kit and the Omni Plant RNA Kit (tDNase I) (Tiangen, Beijing, China)[37].

    cDNA was synthesized for qPCR by using the PrimeScript First Strand cDNA Synthesis Kit (Takara, Dalian, China). The cDNA for qPCR was synthesized by using the PrimeScript First Strand cDNA Synthesis Kit (Takara, Dalian, China). Quantitative real-time fluorescence analysis was performed by using the UltraSYBR Mixture (Low Rox) kit (ComWin Biotech Co. Ltd., Beijing, China). qRT-PCR experiments were performed using the 2−ΔΔCᴛ method for data analysis. The data were analyzed by the 2−ΔΔCᴛ method[31].

    GUS staining buffer contained 1 mM 5-bromo-4-chloro-3-indolyl-β-glutamic acid, 0.01 mM EDTA, 0.5 mM hydrogen ferrocyanide, 100 mM sodium phosphate (pH 7.0), and 0.1% (v/v) Triton X-100 was maintained at 37 °C in the dark. The pMdNRT1.1::GUS construct was transiently introduced into apple calli. To confirm whether MdNRT1.1 is activated or inhibited by MdARF18, we co-transformed 35S::MdARF18 into pMdNRT1.1::GUS is calling. The activity of transgenic calli was assessed using GUS labeling and activity assays[33,38].

    The specimens were crushed into fine particles, combined with 1 mL of ddH2O, and thereafter subjected to a temperature of 100 °C for 30 min. The supernatant was collected in a flow cell after centrifugation at 12,000 revolutions per minute for 10 min. The AutoAnalyzer 3 continuous flow analyzer was utilized to measure nitrate concentrations. (SEAL analytical, Mequon, WI, USA). Nitrate reductase (NR) activity was characterized by the corresponding kits (Solarbio Life Science, Beijing, China) using a spectrophotometric method[31].

    Y1H assays were performed as previously described by Liu et al.[39]. The coding sequence of MdARF18 was integrated into the pGADT7 expression vector, whereas the promoter region of MdNRT1.1 was included in the pHIS2 reporter vector. Subsequently, the constitutive vectors were co-transformed into the yeast monohybrid strain Y187. The individual transformants were assessed on a medium lacking tryptophan, leucine, and histidine (SDT/-L/-H). Subsequently, the positive yeast cells were identified using polymerase chain reaction (PCR). The yeast strain cells were diluted at dilution factors of 10, 100, 1,000, and 10,000. Ten μL of various doses were added to selective medium (SD-T/-L/-H) containing 120 mM 3-aminotriazole (3-AT) and incubated at 28 °C for 2−3 d[37].

    Dual-luciferase assays were performed as described previously[40]. Full-length MdARF18 was cloned into pGreenII 62-SK to produce MdARF18-62-SK. The promoter fragment of MdNRT1.1 was cloned into pGreenII 0800-LUC to produce pMdNRT1.1-LUC. Different combinations were transformed into Agrobacterium tumefaciens LBA4404 and the Agrobacterium solution was injected onto the underside of the leaves of tobacco (Nicotiana benthamiana) leaves abaxially. The Dual Luciferase Reporter Kit (Promega, www.promega.com) was used to detect fluorescence activity.

    Total protein was extracted from wild-type and transgenic apple calli with or without 100 μM MG132 treatment. The purified MdARF18-HIS fusion protein was incubated with total protein[41]. Samples were collected at the indicated times (0, 1, 3, 5, and 7 h).

    Protein gel blots were analyzed using GST antibody. ACTIN antibody was used as an internal reference. All antibodies used in this study were provided by Abmart (www.ab-mart.com).

    Unless otherwise noted, every experiment was carried out independently in triplicate. A one-way analysis of variance (ANOVA) was used to establish the statistical significance of all data, and Duncan's test was used to compare results at the p < 0.05 level[31].

    To investigate whether auxin affects the effective uptake of nitrate in apple, we first externally applied IAA under normal N (5 mM NO3) environment, and this result showed that the growth of Gala apple seedlings in the IAA-treated group were better than the control, and their fresh weights were heavier than the control group (Fig. 1a, d). The N-related physiological indexes of apple seedlings also showed that the nitrate content and NR activity of the root part of the IAA-treated group were significantly higher than the control group, while the nitrate content and NR activity of the shoot part were lower than the control group (Fig. 1b, c). These results demonstrate that auxin could promote the uptake of nitrate and thus promotes growth of plants.

    Figure 1.  Auxin enhances nitrate uptake of Gala seedlings. (a) Phenotypes of apple (Malus domestica 'Royal Gala') seedlings grown nutritionally for 50 d under IAA (10 μM) treatment. (b) Nitrate content of shoot and root apple (Malus domestica 'Royal Gala') seedlings treated with IAA. (c) NR activity in shoot and root of IAA treatment apple (Malus domestica 'Royal Gala') seedlings. (d) Seedling fresh weight under IAA treatment. Bars represent the mean ± SD (n = 3). Different letters above the bars indicate significant differences using the LSD test (p < 0.05).

    To test whether auxin affects the expression of genes related to nitrogen uptake and metabolism. For the root, the expression levels of MdNRT1.1, MdNRT2.1, MdNIA1, MdNIA2, and MdNIR were higher than control group (Supplemental Fig. S1a, f, hj), while the expression levels of MdNRT1.2, MdNRT1.6 and MdNRT2.5 were lower than control group significantly (Supplemental Fig. S1b, d, g). For the shoot, the expression of MdNRT1.1, MdNRT1.5, MdNRT1.6, MdNRT1.7, MdNRT2.1, MdNRT2.5, MdNIA1, MdNIA2, and MdNIR genes were significantly down-regulated (Supplemental Fig. S1a, cj). This result infers that the application of auxin could mediate nitrate uptake in plants by affecting the expression levels of relevant nitrate uptake and assimilation genes.

    Since the auxin signaling pathway requires the regulation of the auxin response factors (ARFs)[25,27], it was investigated whether members of ARF genes were nitrate responsive. Firstly, qPCR quantitative analysis showed that the five subfamily genes of MdARFs (MdARF9, MdARF2, MdARF12, MdARF3, and MdARF18) were expressed at different levels in various organs of the plant (Supplemental Fig. S2). Afterward, the expression levels of five ARF genes were analyzed under different concentrations of nitrate treatment (Fig. 2), and it was concluded that these genes represented by each subfamily responded in different degrees, but the expression level of MdARF18 was up-regulated regardless of low or high nitrogen (Fig. 2i, j), and the expression level of MdARF18 showed a trend of stable up-regulation under IAA treatment (Supplemental Fig. S3). The result demonstrates that MdARFs could affect the uptake of external nitrate by plants and MdARF18 may play an important role in the regulation of nitrate uptake.

    Figure 2.  Relative expression analysis of MdARFs subfamilies in response to different concentrations of nitrate. Expression analysis of representative genes from five subfamilies of MdARF transcription factors. Bars represent the mean ± SD (n = 3). Different letters above the bars indicate significant differences using the LSD test (p < 0.05).

    MdARF18 (MD07G1152100) was predicted through The Apple Genome website (https://iris.angers.inra.fr/gddh13/) and it had high fitness with AtARF18 (AT3G61830). The homologs of ARF18 from 15 species were then identified in NCBI (www.ncbi.nlm.nih.gov) and then constructed an evolutionary tree (Supplemental Fig. S4). The data indicates that MdARF18 was most closely genetically related to MbARF18 (Malus baccata), indicating that they diverged recently in evolution (Supplemental Fig. S4). Conserved structural domain analyses indicated that all 15 ARF18 proteins had highly similar conserved structural domains (Supplemental Fig. S5). In addition, multiple sequence alignment analysis showed that all 15 ARF18 genes have B3-type DNA-binding domains (Supplemental Fig. S6), which is in accordance with the previous reports on ARF18 protein structure[26].

    To explore whether MdARF18 could affect the development of the plant's root system. Firstly, MdARF18 was heterologously expressed into Arabidopsis, and an arf18 mutant (GABI_699B09) Arabidopsis was also obtained (Supplemental Fig. S7). Seven-day-old MdARF18 transgenic Arabidopsis and arf18 mutants were treated in a medium with different nitrate concentrations for 10 d (Fig. 3a, b). After observing results, it was found that under the environment of high nitrate concentration, the primary root of MdARF18 was shorter than arf18 and wild type (Fig. 3c), and the primary root length of arf18 is the longest (Fig. 3c), while there was no significant difference in the lateral root (Fig. 3d). For low nitrate concentration, there was no significant difference in the length of the primary root, and the number of lateral roots of MdARF18 was slightly more than wild type and arf18 mutant. These results suggest that MdARF18 affects root development in plants. However, in general, low nitrate concentrations could promote the transport of IAA by NRT1.1 and thus inhibit lateral root production[3], so it might be hypothesized that MdARF18 would have some effect on MdNRT1.1 thus leading to the disruption of lateral root development.

    Figure 3.  MdARF18 inhibits root development. (a) MdARF18 inhibits root length at 10 mM nitrate concentration. (b) MdARF18 promotes lateral root growth at 0.5 mM nitrate concentration. (c) Primary root length statistics. (d) Lateral root number statistics. Bars represent the mean ± SD (n = 3). Different letters above the bars indicate significant differences using the LSD test (p < 0.05).

    To investigate whether MdARF18 affects the growth of individual plants under different concentrations of nitrate, 7-day-old overexpression MdARF18, and arf18 mutants were planted in the soil and incubated for 20 d. It was found that arf18 had the best growth of shoot, while MdARF18 had the weakest shoot growth at any nitrate concentration (Fig. 4a). MdARF18 had the lightest fresh weight and the arf18 mutant had the heaviest fresh weight (Fig. 4b). N-related physiological indexes revealed that the nitrate content and NR activity of arf18 were significantly higher than wild type, whereas MdARF18 materials were lower than wild type (Fig. 4c, d). More detail, MdARF18 had the lightest fresh weight under low and normal nitrate, while the arf18 mutant had the heaviest fresh weight, and the fresh weight of arf18 under high nitrate concentration did not differ much from the wild type (Fig. 4b). Nitrogen-related physiological indexes showed that the nitrate content of arf18 was significantly higher than wild type, while MdARF18 was lower than wild type. The NR activity of arf18 under high nitrate did not differ much from the wild type, but the NR activity of MdARF18 was the lowest in any treatment (Fig. 4c, d). These results indicate that MdARF18 significantly inhibits plant growth by inhibiting plants to absorb nitrate, and is particularly pronounced at high nitrate concentrations.

    Figure 4.  Ectopic expression of MdARF18 inhibits Arabidopsis growth. (a) Status of Arabidopsis growth after one month of incubation at different nitrate concentrations. (b) Fresh weight of Arabidopsis. (c) Nitrate content of Arabidopsis. (d) NR activity in Arabidopsis. Bars represent the mean ± SD (n = 3). Different letters above the bars indicate significant differences using the LSD test (p < 0.05).

    In addition, to further validate this conclusion, MdARF18 overexpression calli were obtained and treated with different concentrations of nitrate (Supplemental Fig. S8). The results show that the growth of overexpressed MdARF18 was weaker than wild type in both treatments (Supplemental Fig. S9a). The fresh weight of MdARF18 was significantly lighter than wild type (Supplemental Fig. S9b), and its nitrate and NR activity were lower than wild type (Supplemental Fig. S9c, d), which was consistent with the above results (Fig. 4). This result further confirms that MdARF18 could inhibit the development of individual plants by inhibiting the uptake of nitrate.

    Nitrate acts as a signaling molecule that takes up nitrate by activating the NRT family as well as NIAs and NIR[3,34]. To further investigate the pathway by which MdARF18 inhibits plant growth and reduces nitrate content, qRT-PCR was performed on the above plant materials treated with different concentrations of nitrate (Fig. 5). The result shows that the expression levels of AtNRT1.1, AtNIA1, AtNIA2, and AtNIR were all down-regulated in overexpression of MdARF18, and up-regulated in the arf18 mutant (Fig. 5a, hj). There was no significant change in AtNRT1.2 at normal nitrate levels, but AtNRT1.2 expression levels were down-regulated in MdARF18 and up-regulated in arf18 at both high and low nitrate levels (Fig. 5b). This trend in the expression levels of these genes might be consistent with the fact that MdARF18 inhibits the expression of nitrogen-related genes and restricts plant growth. The trend in the expression levels of these genes is consistent with MdARF18 restricting plant growth by inhibiting the expression of nitrogen-related genes. However, AtNRT1.5, AtNRT1.6, AtNRT1.7, AtNRT2.1, and AtNRT2.5 did not show suppressed expression levels in MdARF18 (Fig. 5cg). These results suggest that MdARF18 inhibits nitrate uptake and plant growth by repressing some of the genes for nitrate uptake or assimilation.

    Figure 5.  qPCR-RT analysis of N-related genes. Expression analysis of N-related genes in MdARF18 transgenic Arabidopsis at different nitrate concentrations. Bars represent the mean ± SD (n = 3). Different letters above the bars indicate significant differences using the LSD test (p < 0.05).

    In addition, to test whether different concentrations of nitrate affect the protein stability of MdARF18. However, it was found that there was no significant difference in the protein stability of MdARF18 at different concentrations of nitrate (Supplemental Fig. S10). This result suggests that nitrate does not affect the degradation of MdARF18 protein.

    To further verify whether MdARF18 can directly bind N-related genes, firstly we found that the MdNRT1.1 promoter contains binding sites to ARF factors (Fig. 6a). The yeast one-hybrid research demonstrated an interaction between MdARF18 and the MdNRT1.1 promoter, as shown in Fig. 6b. Yeast cells that were simultaneously transformed with MdNRT1.1-P-pHIS and pGADT7 were unable to grow in selected SD medium. However, cells that were transformed with MdNRT1.1-P-pHIS and MdARF18-pGADT7 grew successfully in the selective medium. The result therefore hypothesizes that MdARF18 could bind specifically to MdNRT1.1 promoter to regulate nitrate uptake in plants.

    Figure 6.  MdARF18 binds directly to the promoter of MdNRT1.1. (a) Schematic representation of MdNRT1.1 promoter. (b) Y1H assay of MdARF18 bound to the MdNRT1.1 promoter in vitro. 10−1, 10−2, 10−3, and 10−4 indicate that the yeast concentration was diluted 10, 100, 1,000, and 10,000 times, respectively. 3-AT stands for 3-Amino-1,2,4-triazole. (c) Dual luciferase assays demonstrate the binding of MdARF18 with MdNRT1.1 promoter. The horizontal bar on the left side of the right indicates the captured signal intensity. Empty LUC and 35S vectors were used as controls. Representative images of three independent experiments are shown here.

    To identify the inhibition or activation of MdNRT1.1 by MdARF18, we analyzed their connections by Dual luciferase assays (Fig. 6c), and also analyzed the fluorescence intensity (Supplemental Fig. S11). It was concluded that the fluorescence signals of cells carrying 35Spro and MdNRT1.1pro::LUC were stronger, but the mixture of 35Spro::MdARF18 and MdNRT1.1pro::LUC injected with fluorescence signal intensity was significantly weakened. Next, we transiently transformed the 35S::MdARF18 into pMdNRT1.1::GUS transgenic calli (Fig. 7). GUS results first showed that the color depth of pMdNRT1.1::GUS and 35S::MdARF18 were significantly lighter than pMdnNRT1.1::GUS alone (Fig. 7a). GUS enzyme activity, as well as GUS expression, also indicated that the calli containing pMdnNRT1.1::GUS alone had a stronger GUS activity (Fig. 7b, c). In addition, the GUS activity of calli containing both pMdNRT1.1:GUS and 35S::MdARF18 were further attenuated under both high and low nitrate concentrations (Fig. 7a). These results suggest that MdARF18 represses MdNRT1.1 expression by directly binding to the MdNRT1.1 promoter region.

    Figure 7.  MdARF18 inhibits the expression of MdNRT1.1. (a) GUS staining experiment of pMdNRT1.1::GUS transgenic calli and transgenic calli containing both pMdNRT1.1::GUS and 35S::MdARF18 with different nitrate treatments. (b) GUS activity assays in MdARF18 overexpressing calli with different nitrate treatments. (c) GUS expression level in MdARF18 overexpressing calli with different nitrate treatments. Bars represent the mean ± SD (n = 3). Different numbers of asterisk above the bars indicate significant differences using the LSD test (*p < 0.05 and **p< 0.01).

    Plants replenish their nutrients by absorbing nitrates from the soil[42,43]. Previous studies have shown that some of the plant hormones such as IAA, GA, and ABA interact with nitrate[25,4445]. The effect of nitrate on the content and transport of IAA has been reported in previous studies, e.g., nitrate supply reduced IAA content in Arabidopsis, wheat, and maize roots and inhibited the transport of IAA from shoot to root[20,21]. In this study, it was found that auxin treatment promoted individual fresh weight gain and growth (Fig. 1a, b). Nitrate content and NR activity were also significantly higher in their root parts (Fig. 1c, d) and also affected the transcript expression levels of related nitrate uptake and assimilation genes (Supplemental Fig. S1). Possibly because IAA can affect plant growth by influencing the uptake of external nitrates by the plant.

    ARFs are key transcription factors to regulate auxin signaling[4649]. We identified five representative genes of the apple MdARFs subfamily and they all had different expression patterns (Supplemental Fig. S2). The transcript levels of each gene were found to be affected to different degrees under different concentrations of nitrate, but the expression level of MdARF18 was up-regulated under both low and high nitrate conditions (Fig. 2). The transcript level of MdARF18 was also activated under IAA treatment (Supplemental Fig. S3), so MdARF18 began to be used in the study of the mechanism of nitrate uptake in plants. In this study, an Arabidopsis AtARF18 homolog was successfully cloned and named MdARF18 (Supplemental Figs S4, S5). It contains a B3-type DNA-binding structural domain consistent with previous studies of ARFs (Supplemental Fig. S6), and arf18 mutants were also obtained and their transcript levels were examined (Supplemental Fig. S7).

    Plants rely on rapid modification of the root system to efficiently access effective nitrogen resources in the soil for growth and survival. The plasticity of root development is an effective strategy for accessing nitrate, and appropriate concentrations of IAA can promote the development of lateral roots[7,44]. The present study found that the length of the primary root was shortened and the number of lateral roots did increase in IAA-treated Gl3 apple seedlings (Supplemental Fig. S12). Generally, an environment with low concentrations of nitrate promotes the transport of IAA by AtNRT1.1, which inhibits the growth of lateral roots[14]. However, in the research of MdARF18 transgenic Arabidopsis, it was found that the lateral roots of MdARF18-OX increased under low concentrations of nitrate, but there was no significant change in the mutant arf18 (Fig. 3). Therefore, it was hypothesized that MdARF18 might repress the expression of the MdNRT1.1 gene or other related genes that can regulate root plasticity, thereby affecting nitrate uptake in plants.

    In rice, several researchers have demonstrated that OsARF18 significantly regulates nitrogen utilization. Loss of function of the Rice Salt Tolerant 1 (RST1) gene (encoding OsARF18) removes its ability to transcriptionally repress OsAS1, accelerating the assimilation of NH4+ to Asn and thereby increasing nitrogen utilization[28]. During soil incubation of MdARF18-OX Arabidopsis, it was found that leaving aside the effect of differences in nitrate concentration, the arf18 mutant grew significantly better than MdARF18-OX and had higher levels of nitrate and NR activity in arf18 than in MdARF18-OX. This demonstrates that MdARF18 may act as a repressor of nitrate uptake and assimilation, thereby inhibiting normal plant development (Fig. 4). Interestingly, an adequate nitrogen environment promotes plant growth, but MdARF18-OX Arabidopsis growth and all physiological indexes were poorer under high nitrate concentration than MdARF18-OX at other concentrations. We hypothesize that MdARF18 may be activated more intensively at high nitrate concentrations, or that MdARF18 suppresses the expression levels of genes for nitrate uptake or assimilation (genes that may play a stronger role at high nitrate concentrations), thereby inhibiting plant growth. In addition, we obtained MdARF18 transgenic calli (Supplemental Fig. S8) and subjected them to high and low concentrations of nitrate, and also found that MdARF18 inhibited the growth of individuals at both concentrations (Supplemental Fig. S9). This further confirms that MdARF18 inhibits nitrate uptake in individuals.

    ARF family transcription factors play a key role in transmitting auxin signals to alter plant growth and development, e.g. osarf1 and osarf24 mutants have reduced levels of OsNRT1.1B, OsNRT2.3a and OsNIA2 transcripts[22]. Therefore, further studies are needed to determine whether MdARF18 activates nitrate uptake through different molecular mechanisms. The result revealed that the transcript levels of AtNRT1.1, AtNIA1, AtNIA2, and AtNIR in MdARF18-OX were consistent with the developmental pattern of impaired plant growth (Fig. 5). Unfortunately, we attempted to explore whether variability in nitrate concentration affects MdARF18 to differ at the protein level, but the two did not appear to differ significantly (Supplemental Fig. S10).

    ARF transcription factors act as trans-activators/repressors of N metabolism-related genes by directly binding to TGTCNN/NNGACA-containing fragments in the promoter regions of downstream target genes[27,50]. The NRT family plays important roles in nitrate uptake, transport, and storage, and NRT1.1 is an important dual-affinity nitrate transporter protein[7,5052], and nitrogen utilization is very important for apple growth[53,54]. We identified binding sites in the promoters of these N-related genes that are compatible with ARF factors, and MdARF18 was found to bind to MdNRT1.1 promoter by yeast one-hybrid technique (Fig. 6a, b). It was also verified by Dual luciferase assays that MdARF18 could act as a transcriptional repressor that inhibited the expression of the downstream gene MdNRT1.1 (Fig. 6c), which inhibited the uptake of nitrate in plants. In addition, the GUS assay was synchronized to verify that transiently expressed pMdNRT1.1::GUS calli with 35S::MdARF18 showed a lighter staining depth and a significant decrease in GUS transcript level and enzyme activity (Fig. 7). This phenomenon was particularly pronounced at high concentrations of nitrate. These results suggest that MdARF18 may directly bind to the MdNRT1.1 promoter and inhibit its expression, thereby suppressing NO3 metabolism and decreasing the efficiency of nitrate uptake more significantly under high nitrate concentrations.

    In conclusion, in this study, we found that MdARF18 responds to nitrate and could directly bind to the TGTCTT site of the MdNRT1.1 promoter to repress its expression. Our findings provide new insights into the molecular mechanisms by which MdARF18 regulates nitrate transport in apple.

    The authors confirm contribution to the paper as follows: study conception and design: Liu GD; data collection: Liu GD, Rui L, Liu RX; analysis and interpretation of results: Liu GD, Li HL, An XH; draft manuscript preparation: Liu GD; supervision: Zhang S, Zhang ZL; funding acquisition: You CX, Wang XF; All authors reviewed the results and approved the final version of the manuscript.

    Data sharing not applicable to this article as no datasets were generated or analyzed during the current study.

    This work was supported by the National Natural Science Foundation of China (32272683), the Shandong Province Key R&D Program of China (2022TZXD008-02), the China Agriculture Research System of MOF and MARA (CARS-27), the National Key Research and Development Program of China (2022YFD1201700), and the National Natural Science Foundation of China (NSFC) (32172538).

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

  • Supplemental File 1 List of differentially expressed genes (DEGs) in different comparisons.
    Supplemental File 2 Gene Ontology (GO) classification of the DEGs identified in different comparisons.
    Supplemental File 3 KEGG classification of the DEGs identified in different comparisons.
    Supplemental File 4 List of candidate acclimation response genes that be identified both in QTL regions and differentially expressed genes. QTL regions were previously reported in Brown et al.[15].
    Supplemental File 5 List of candidate acclimation response genes that were also differentially expressed in protein abundance in a previous study[16].
  • [1]

    Li M, Yuyama N, Hirata M, Han J, Wang Y, et al. 2009. Construction of a high-density SSR marker-based linkage map of zoysiagrass (Zoysia japonica Steud.). Euphytica 170:327−38

    doi: 10.1007/s10681-009-9990-8

    CrossRef   Google Scholar

    [2]

    Patton AJ, Schwartz BM, Kenworthy KE. 2017. Zoysiagrass (Zoysia spp.) history, utilization, and improvement in the United States: a review. Crop Science 57:S-37−S-72

    doi: 10.2135/cropsci2017.02.0074

    CrossRef   Google Scholar

    [3]

    Hinton JD, Livingston DP, III, Miller GL, Peacock CH, Tuong T. 2012. Freeze tolerance of nine zoysiagrass cultivars using natural cold acclimation and freeze chambers. HortScience 47:112−15

    doi: 10.21273/HORTSCI.47.1.112

    CrossRef   Google Scholar

    [4]

    Patton AJ, Cunningham SM, Volenec JJ, Reicher ZJ. 2007. Differences in freeze tolerance of zoysiagrasses: II. Carbohydrate and proline accumulation. Crop Science 47:2170−81

    doi: 10.2135/cropsci2006.12.0784

    CrossRef   Google Scholar

    [5]

    Baier M, Bittner A, Prescher A, van Buer J. 2019. Preparing plants for improved cold tolerance by priming. Plant, Cell & Environment 42:782−800

    doi: 10.1111/pce.13394

    CrossRef   Google Scholar

    [6]

    Chinnusamy V, Zhu J, Zhu J. 2007. Cold stress regulation of gene expression in plants. Trends in Plant Science 12:444−51

    doi: 10.1016/j.tplants.2007.07.002

    CrossRef   Google Scholar

    [7]

    Huang X, Cao L, Fan J, Ma G, Chen L. 2022. CdWRKY2-mediated sucrose biosynthesis and CBF-signalling pathways coordinately contribute to cold tolerance in bermudagrass. Plant Biotechnology Journal 20:660−75

    doi: 10.1111/pbi.13745

    CrossRef   Google Scholar

    [8]

    Gilmour SJ, Sebolt AM, Salazar MP, Everard JD, Thomashow MF. 2000. Overexpression of the Arabidopsis CBF3 transcriptional activator mimics multiple biochemical changes associated with cold acclimation. Plant Physiology 124:1854−65

    doi: 10.1104/pp.124.4.1854

    CrossRef   Google Scholar

    [9]

    Gilmour SJ, Fowler SG, Thomashow MF. 2004. Arabidopsis transcriptional activators CBF1, CBF2, and CBF3 have matching functional activities. Plant Molecular Biology 54:767−81

    doi: 10.1023/B:PLAN.0000040902.06881.d4

    CrossRef   Google Scholar

    [10]

    Wei S, Du Z, Gao F, Ke X, Li J, et al. 2015. Global transcriptome profiles of 'Meyer' Zoysiagrass in response to cold stress. PLoS ONE 10:e0131153

    doi: 10.1371/journal.pone.0131153

    CrossRef   Google Scholar

    [11]

    Chinnusamy V, Ohta M, Kanrar S, Lee BH, Hong X, et al. 2003. ICE1: a regulator of cold-induced transcriptome and freezing tolerance in Arabidopsis. Genes & Development 17:1043−54

    doi: 10.1101/gad.1077503

    CrossRef   Google Scholar

    [12]

    Hundertmark M, Hincha DK. 2008. LEA (Late Embryogenesis Abundant) proteins and their encoding genes in Arabidopsis thaliana. BMC Genomics 9:118

    doi: 10.1186/1471-2164-9-118

    CrossRef   Google Scholar

    [13]

    Thomashow MF. 1998. Role of cold-responsive genes in plant freezing tolerance. Plant Physiology 118:1−8

    doi: 10.1104/pp.118.1.1

    CrossRef   Google Scholar

    [14]

    Holloway HMP, Yu X, Dunne JC, Schwartz BM, Patton AJ, et al. 2018. A SNP-based high-density linkage map of zoysiagrass (Zoysia japonica Steud.) and its use for the identification of QTL associated with winter hardiness. Molecular Breeding 38:10

    doi: 10.1007/s11032-017-0763-0

    CrossRef   Google Scholar

    [15]

    Brown JM, Yu X, Holloway HMP, Tuong TD, Schwartz BM, et al. 2021. Identification of QTL associated with cold acclimation and freezing tolerance in Zoysia japonica. Crop Science 61:3044−55

    doi: 10.1002/csc2.20368

    CrossRef   Google Scholar

    [16]

    Brown JM, Yu X, Holloway HMP, DaCosta M, Bernstein RP, et al. 2020. Differences in proteome response to cold acclimation in Zoysia japonica cultivars with different levels of freeze tolerance. Crop Science 60:2744−56

    doi: 10.1002/csc2.20225

    CrossRef   Google Scholar

    [17]

    Han Q, Zhu Q, Shen Y, Lee M, Lübberstedt T, et al. 2022. QTL mapping low-temperature germination ability in the maize IBM Syn10 DH population. Plants 11:214

    doi: 10.3390/plants11020214

    CrossRef   Google Scholar

    [18]

    Lei L, Zheng H, Bi Y, Yang L, Liu H, et al. 2020. Identification of a major QTL and candidate gene analysis of salt tolerance at the bud burst stage in rice (Oryza sativa L.) using QTL-Seq and RNA-Seq. Rice 13:55

    doi: 10.1186/s12284-020-00416-1

    CrossRef   Google Scholar

    [19]

    Patton AJ. 2009. Selecting zoysiagrass cultivars: turfgrass quality, growth, pest and environmental stress tolerance. Applied Turfgrass Science 6:1−18

    doi: doi.org/10.1094/ats-2009-1019-01-mg

    CrossRef   Google Scholar

    [20]

    Chen L, Fan J, Hu L, Hu Z, Xie Y, et al. 2015. A transcriptomic analysis of bermudagrass (Cynodon dactylon) provides novel insights into the basis of low temperature tolerance. BMC Plant Biology 15:216

    doi: 10.1186/s12870-015-0598-y

    CrossRef   Google Scholar

    [21]

    Anderson JA, Taliaferro CM, Martin DL. 1993. Evaluating freeze tolerance of bermudagrass in a controlled environment. HortScience 28:955

    doi: 10.21273/HORTSCI.28.9.955

    CrossRef   Google Scholar

    [22]

    Tanaka H, Hirakawa H, Kosugi S, Nakayama S, Ono A, et al. 2016. Sequencing and comparative analyses of the genomes of zoysiagrasses. DNA Research 23:171−80

    doi: 10.1093/dnares/dsw006

    CrossRef   Google Scholar

    [23]

    Anders S, Pyl PT, Huber W. 2015. HTSeq—a Python framework to work with high-throughput sequencing data. Bioinformatics 31:166−69

    doi: 10.1093/bioinformatics/btu638

    CrossRef   Google Scholar

    [24]

    Anders S, Huber W. 2010. Differential expression analysis for sequence count data. Genome Biology 11:R106

    doi: 10.1186/gb-2010-11-10-r106

    CrossRef   Google Scholar

    [25]

    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

    [26]

    Holloway HPM. 2016. Genomic approaches for improving freeze tolerance in zoysiagrass. Thesis. North Carolina State University, U. S.

    [27]

    Guy CL. 1990. Cold acclimation and freezing stress tolerance: role of protein metabolism. Annual Review of Plant Physiology and Plant Molecular Biology 41:187−223

    doi: 10.1146/annurev.pp.41.060190.001155

    CrossRef   Google Scholar

    [28]

    Hannah MA, Heyer AG, Hincha DK. 2005. A global survey of gene regulation during cold acclimation in Arabidopsis thaliana. PLoS Genetics 1:e26

    doi: 10.1371/journal.pgen.0010026

    CrossRef   Google Scholar

    [29]

    Maruyama K, Sakuma Y, Kasuga M, Ito Y, Seki M, et al. 2004. Identification of cold-inducible downstream genes of the Arabidopsis DREB1A/CBF3 transcriptional factor using two microarray systems. The Plant Journal 38:982−93

    doi: 10.1111/j.1365-313X.2004.02100.x

    CrossRef   Google Scholar

    [30]

    Liu Z, Jia Y, Ding Y, Shi Y, Li Z, et al. 2017. Plasma membrane CRPK1-mediated phosphorylation of 14-3-3 proteins induces their nuclear import to fine-tune CBF signaling during cold response. Molecular Cell 66:117−128.E5

    doi: 10.1016/j.molcel.2017.02.016

    CrossRef   Google Scholar

    [31]

    Finkel T. 2011. Signal transduction by reactive oxygen species. Journal of Cell Biology 194:7−15

    doi: 10.1083/jcb.201102095

    CrossRef   Google Scholar

    [32]

    Lin F, Zhou L, He B, Zhang X, Dai H, et al. 2019. QTL mapping for maize starch content and candidate gene prediction combined with co-expression network analysis. Theoretical and Applied Genetics 132:1931−41

    doi: 10.1007/s00122-019-03326-z

    CrossRef   Google Scholar

    [33]

    Septiani P, Lanubile A, Stagnati L, Busconi M, Nelissen H, et al. 2019. Unravelling the genetic basis of Fusarium seedling rot resistance in the MAGIC maize population: novel targets for breeding. Scientific Reports 9:5665

    doi: 10.1038/s41598-019-42248-0

    CrossRef   Google Scholar

    [34]

    Naoumkina M, Thyssen GN, Fang DD, Jenkins JN, Mccarty JC, et al. 2019. Genetic and transcriptomic dissection of the fiber length trait from a cotton (Gossypium hirsutum L.) MAGIC population. BMC Genomics 20:112

    doi: 10.1186/s12864-019-5427-5

    CrossRef   Google Scholar

    [35]

    Kim KN, Cheong YH, Gupta R, Luan S. 2000. Interaction specificity of Arabidopsis calcineurin B-Like calcium sensors and their target kinases. Plant Physiology 124:1844−53

    doi: 10.1104/pp.124.4.1844

    CrossRef   Google Scholar

    [36]

    Huang C, Ding S, Zhang H, Du H, An L. 2011. CIPK7 is involved in cold response by interacting with CBL1 in Arabidopsis thaliana. Plant Science 181:57−64

    doi: 10.1016/j.plantsci.2011.03.011

    CrossRef   Google Scholar

    [37]

    Wu C. 1995. Heat shock transcription factors: structure and regulation. Annual Review of Cell and Developmental Biology 11:441−69

    doi: 10.1146/annurev.cb.11.110195.002301

    CrossRef   Google Scholar

    [38]

    Sasaki K, Christov NK, Tsuda S, Imai R. 2014. Identification of a novel LEA protein involved in freezing tolerance in wheat. Plant and Cell Physiology 55:136−47

    doi: 10.1093/pcp/pct164

    CrossRef   Google Scholar

    [39]

    Ogaya R, Peñuelas J, Asensio D, Llusià J. 2011. Chlorophyll fluorescence responses to temperature and water availability in two co-dominant Mediterranean shrub and tree species in a long-term field experiment simulating climate change. Environmental and Experimental Botany 71:123−27

    doi: 10.1016/j.envexpbot.2010.10.016

    CrossRef   Google Scholar

    [40]

    Zhang J, Zhang Z, Liu W, Li L, Han L, et al. 2022. Transcriptome analysis revealed a positive role of ethephon on Chlorophyll metabolism of Zoysia japonica under cold stress. Plants 11:442

    doi: 10.3390/plants11030442

    CrossRef   Google Scholar

    [41]

    Li H, Ding Y, Shi Y, Zhang X, Zhang S, et al. 2017. MPK3- and MPK6-mediated ICE1 phosphorylation negatively regulates ICE1 stability and freezing tolerance in Arabidopsis. Developmental Cell 43:630−642.E4

    doi: h10.1016/j.devcel.2017.09.025

    CrossRef   Google Scholar

    [42]

    Zhao C, Wang P, Si T, Hsu CC, Wang L, et al. 2017. MAP kinase cascades regulate the cold response by modulating ICE1 protein stability. Developmental Cell 43:618−629.E5

    doi: 10.1016/j.devcel.2017.09.024

    CrossRef   Google Scholar

    [43]

    Kong Q, Qu N, Gao M, Zhang Z, Ding X, et al. 2012. The MEKK1-MKK1/MKK2-MPK4 kinase cascade negatively regulates immunity mediated by a mitogen-activated protein kinase kinase kinase in Arabidopsis. The Plant Cell 24:2225−36

    doi: 10.1105/tpc.112.097253

    CrossRef   Google Scholar

  • Cite this article

    Brown JM, Weldt CE, Holloway HMP, Tuong TD, Patton AJ, et al. 2023. Transcriptomic analysis of zoysiagrass (Zoysia japonica) provides novel insights into the molecular basis of cold acclimation. Grass Research 3:25 doi: 10.48130/GR-2023-0025
    Brown JM, Weldt CE, Holloway HMP, Tuong TD, Patton AJ, et al. 2023. Transcriptomic analysis of zoysiagrass (Zoysia japonica) provides novel insights into the molecular basis of cold acclimation. Grass Research 3:25 doi: 10.48130/GR-2023-0025

Figures(5)  /  Tables(3)

Article Metrics

Article views(3982) PDF downloads(591)

ARTICLE   Open Access    

Transcriptomic analysis of zoysiagrass (Zoysia japonica) provides novel insights into the molecular basis of cold acclimation

Grass Research  3 Article number: 25  (2023)  |  Cite this article

Abstract: Zoysia japonica (Z. japonica) is a warm-season perennial turfgrass commonly grown in the southeastern United States for its relatively low input requirements and general tolerance to drought, shade, and salinity. Improving freezing tolerance is critical for Z. japonica, as it could extend the northern boundaries that the species is able to be grown in. To deepen our understanding of the molecular basis of freezing tolerance in Z. japonica, a transcriptomic approach was utilized to identify genes involved in cold acclimation. 'Meyer', a freeze tolerant cultivar, and 'Victoria', a freeze susceptible cultivar were subjected to cold acclimation and non-cold acclimation treatments to determine the number of differentially expressed genes (DEGs). In response to cold acclimation, a total of 4,609 DEGs were up-regulated and 3,605 DEGs were down-regulated in 'Meyer', while 3,758 DEGs were up-regulated and 3,516 DEGs were down-regulated in 'Victoria'. GO and KEGG enrichment analysis showed several diverse pathways and biological processes, including photosynthesis, transmembrane transport, and plant hormone signal transduction. Integrating these information with previous studies on proteomics and QTL mapping, several candidate genes were identified to be associated with cold acclimation and freezing tolerance across different studies, such as LEA, CIPK, POD, HSF, HSP, MPK, PSII and multiple transcription factors. The candidate genes identified in this study showed great value in potentially being a target for future selection efforts for freeze-tolerant Z. japonica.

    • Z. japonica is a warm-season turfgrass desirable for both residential and commercial turf use due to its fine leaf texture, low growth habit, low fertility requirements, and wide general tolerance to abiotic stressors such as drought, shade and salinity[1,2]. Although Z. japonica has the best freeze tolerance among warm-season grasses[2], currently, the commercial spread of Z. japonica is still limited to the transition zone and the southern US. In order to expand the species to northern cooler regions, the increased freeze-tolerant cultivars are critically demanded in Z. japonica breeding. Although several cultivars have been tested for cold acclimation and freeze tolerance[3], the lack of knowledge on the molecular mechanisms responsible for cold and freezing response limits breeding progress. Therefore, elucidating the molecular mechanisms and identifying the candidate genes underlying cold stress will accelerate the speed of improvement of freezing tolerance in Z. japonica.

      Cold acclimation is a critical process in the development of freezing tolerance, and is defined as exposure to low but above freezing temperatures prior to a freeze stress[4]. Multiple biochemical changes have been identified as associated with cold acclimation, including an increased generation of reactive oxygen species (ROS), an up-regulation of antioxidant enzymes, alteration in plant hormone levels, and accumulation of amino acids with osmotic effects[5]. Meanwhile, cold acclimation and cold tolerance are involved in complicated regulatory pathways. Numerous genes have been found up- or down-regulated during cold acclimation treatment. Among several cold signaling pathways, the CBF/DREB1-dependent cold signaling pathway is well characterized as the key regulatory pathway[6,7]. Three CBF/DREB1s are involved in the regulation of cold response gene (COR) expression and cold tolerance in Arabidopsis and zoysiagrass[810]. INDUCER OF CBF EXPRESSION (ICE1), a MYC-type bHLH transcription factor, could activate the expression of CBF genes by directly binding to their promoters under cold stress[8]. LEA proteins are thought to be important for membrane stabilization and to prevent protein aggregation during stress[11]. There are several LEA gene members, COR78/RD29A, COR47, COR15A and COR6.6 in Arabidopsis that are induced by cold stress[12,13]. Although many genes and pathways have been identified in model plant species, the molecular mechanism of cold tolerance in Z. japonica is largely unknown.

      Our previous studies have demonstrated that Z. japonica cultivars 'Meyer' and 'Victoria' as freeze tolerant and freeze susceptible genotypes, separately. Holloway et al. reported that 'Victoria' showed more severe winter injury than 'Meyer' in field evaluations[14]. In a lab-based evaluation, 'Meyer' showed less surviving green tissue and regrowth rate than 'Victoria' after freezing without cold acclimation. Meanwhile, the results were reversed if a cold acclimation treatment was added before freezing[15]. In addition, based on LT50 (50% mortality) evaluation, 'Meyer' showed higher freeze tolerance than 'Victoria' after cold acclimation treatment[16]. Meanwhile, genetic analysis identified multiple quantitative trait loci (QTL) regions associated with winter injury, post-freezing surviving green tissue (SGT) and regrowth (RG) using a 'Meyer x Victoria' pseudo-F2 mapping population[14,15]. Several proteins of interest for their association with cold acclimation were identified through a comparative proteomics analysis in 'Meyer' and 'Victoria'[16].

      In recent years, the integration of RNA-seq and QTL mapping has been considered to be a reliable method to rapidly identify potential candidate genes, which has been applied in several crop species like rice and maize[17,18]. However, the strategy has not been reported in turfgrass species. The objectives of this study were: (1) to identify differentially expressed genes (DEGs) and functional annotation in response to cold acclimation in Z. japonica, and (2) to identify potential candidate genes crossing transcriptome with previous proteome and QTL mapping studies. The identified potential candidate genes will provide insight into the molecular mechanisms of cold tolerance, and would be helpful for Z. japonica breeding programs.

    • Z. japonica cultivar 'Meyer' and 'Victoria', have been extensively reported in the literature as freeze tolerant and freeze susceptible genotypes, respectively[4,1416,19], were used in this study. Two-inch plugs of each genetically uniform cultivar were propagated with Fafard 4P potting mix and grown under greenhouse conditions at day/night 27 ± 2 / 22 ± 2 °C for six weeks until fully established. Plants were mowed at 5 cm and fertilized biweekly with 24:8:16 NPK solution of Miracle-Gro Water Soluble All Purpose Plant Food (The Scotts Miracle-Gro® Company, Marysville, OH, USA).

    • The protocol of the cold acclimation treatment was applied according to Chen et al. with minor modifications[20]. Following establishment in the greenhouse, pots were moved into temperature-controlled Percival plant growth chambers (model I-36LL) for cold acclimation. The growth chambers contained a consistent light emitting diode (LED) which maintained a 10-h photoperiod of 300 µmol·m−2·s−1 active radiation[21]. Plants in the cold acclimation (CA) group were held under 8/2 °C day and night temperatures and 12/12 h light cycles for 48 h, while plants in the non-acclimation (NA) group were maintained at a consistent 22 °C with 12/12 h day and night light cycles. After cold acclimation/non-acclimation, 12 leaf samples of similar maturity status were collected for each 'Meyer' and 'Victoria' with three biological replicates. Leaf samples were frozen in liquid nitrogen immediately and stored in a −80 °C freezer until RNA extraction. The experiment contains four treatments, non-cold acclimated 'Meyer' (NM), cold acclimated 'Meyer' (CM), non-cold acclimated 'Victoria' (NV), and cold acclimated 'Victoria' (CV).

    • Total RNA was extracted from frozen leaf samples using the RNeasy Plant Mini Kit (Qiagen, Germany). Only RNA samples that passed quality and quantity checks (RIN value ≥ 7; total amount ≥ 1.5 ug) were used in sequencing library construction. Following this quality control, mRNA was enriched via oligo (dT) beads, and fragmented in fragmentation buffer. First strand cDNA was synthesized using random hexamers and reverse transcriptase, then second strand synthesis buffer (Illumina, San Diego, CA, USA) was added along with dNTPs, RNase H and Escherichia coli polymerase I for second strand synthesis. After a round of purification, terminal repair and poly-A tailing, Illumina PE adapters were ligated and 200 bp cDNA fragments were preferentially selected. Then, cDNA fragments with ligated adapters were enriched using PCR primers. Library quality was assessed by quantification with a Qubit 2.0 fluorometer and qPCR, and insert size was checked on an Agilent 2100. The library was sequenced using an Illumina NovaSeq sequencing platform to generate 150 bp pair-end reads. All the raw sequencing data were deposited into the database of sequence read archive (SRA) in National Center for Biotechnology Information (NCBI) under BioProject ID: PRJNA1040203.

    • The raw reads were processed to remove reads containing adapters, uncertain nucleotides (N > 10%), and low-quality reads. The resulting clean reads were evaluated on Q20, Q30, GC content, and used for downstream analyses. Clean reads were aligned to the Zoysia japonica 'Nagirizaki' transcript using the HISAT2[22]. The FPKM (expected number of Fragments Per Kilo bases of transcript sequence per Millions mapped reads) method was used to estimate gene expression levels by counting reads that mapped to genes or exons using HTSeq v0.6.1[23]. Differential gene expression analysis was performed using the DESeq R package[24], and the significant differential expression was filtered using padj ≤ 0.05 and | log2(Fold Change) | ≥ 1.

    • Differentially expressed genes (DEGs) were assigned to Gene Ontology (GO) enrichment analysis via the GOseq R package[25]. GO with corrected p-values < 0.05 are significantly enriched in DEGs. The KOBAS software was used to test the enrichment of DEGs in KEGG pathways, identifying significantly (corrected p-values < 0.05) enriched metabolic or signal transduction pathways associated with the DEGs.

    • Eleven overlapping QTL regions were identified in our previous studies that were associated with winter injury in field evaluations and with laboratory-based freezing test traits[14,15,26]. We name these regions as QTL1-11 following the chromosome order. The sequences of QTL regions were extracted from the Zoysia japonica 'Nagirizaki' reference genome based on the position of flanking markers[22]. Then, the sequences were blasted against the Zoysia japonica 'Nagirizaki' transcript to obtain the gene ID[22]. Identified genes that overlapped with the above DEGs were recognized as candidate acclimation response genes.

    • A cDNA library was constructed from total RNA of Z. japonica leaf samples under different cold acclimation treatments. Following Illumina sequencing, the overview of sequencing reads obtained is presented in Table 1. Averaged from three replicates for each sample, there were 50,275,906, 60,245,453, 57,550,037, and 60,264,472 raw reads generated from the NM, CM, NV, CV samples, respectively. The raw reads were processed to remove reads containing adapters, uncertain nucleotides (N > 10%), and low-quality reads. After the quality control, average of 49,055,750, 58,957,880, 56,263,343, and 59,058,380 clean reads were obtained for the NM, CM, NV, CV samples, respectively (Table 1).

      Table 1.  Summary of reads following Illumina sequencing.

      SampleRaw readsClean readsClean basesError rate (%)Q20 (%)Q30 (%)GC content (%)
      NM149,475,25448,077,3827.2G0.0397.9294.1557.25
      NM251,183,83650,010,4727.5G0.0397.8493.9957.49
      NM350,168,62849,079,3967.4G0.0397.8894.0457.49
      CM159,804,05658,665,4468.8G0.0397.3792.9552.47
      CM267,204,25865,599,3809.8G0.0397.3592.9451.48
      CM353,728,04652,608,8147.9G0.0397.1192.4852.52
      NV171,039,90269,474,29610.4G0.0397.4793.1953.68
      NV246,536,82445,397,7126.8G0.0397.3993.0352.99
      NV355,073,38453,918,0228.1G0.0397.0892.4153.03
      CV161,931,4586,0855,1049.1G0.0397.4393.1252.9
      CV251,018,63449,819,7847.5G0.0397.2892.8452.43
      CV367,843,32466,500,25210.0G0.0397.4993.2352.78
      *NM: 'Meyer' genotype under non-cold acclimated conditions; CM: 'Meyer' genotype under cold acclimated conditions; NV: 'Victoria' genotype under non-cold acclimated conditions; CV: 'Victoria' genotype under cold acclimated conditions.
    • Mapping results showed that 92.1% of NM, 85.7% of CM, 85.2% of NV and 83.8% of CV clean reads could be mapped to reference genome, which result in a total of 27,209 genes were identified. We paid particular attention to the DEGs identified for comparison of with/without cold acclimation treatment. The comparisons included non-cold acclimated 'Meyer' vs cold acclimated 'Meyer' (R genotype) (NM vs CM) and non-cold acclimated 'Victoria' vs cold acclimated 'Victoria' (S genotype) (NV vs CV). For NM vs CM, 4,609 DEGs were up-regulated and 3,605 DEGs were down-regulated after treatment. For NV vs CV, 3,758 DEGs were up-regulated and 3,516 DEGs were down-regulated in response to cold acclimation (Fig. 1a, Supplemental File 1).

      Figure 1. 

      Differentially expressed genes (DEGs) in Z. japonica in response to cold acclimation. (a) Number of DEGs identified in the 'Meyer' (resistant) and 'Victoria' (susceptible) genotypes. (b) Venn diagram of differentially expressed genes, showing the unique and overlapping genes expressed among comparison groups. NM vs CM: non-cold acclimated 'Meyer' vs cold acclimated 'Meyer' (R genotype); NV vs CV: non-cold acclimated 'Victoria' vs cold acclimated 'Victoria' (S genotype). G1−G8 indicated eight groups.

      To further evaluate DEGs patterns, a Venn diagram was used to show the unique and overlapping genes expressed among comparison groups (Fig. 1b). Accordingly, the DEGs were categorized into eight groups (G1-G8): G1, 2,308 DEGs were up-regulated in both R and S genotypes; G2, 1,631 DEGs were down-regulated in both genotypes; G3, 104 DEGs were up-regulated in R but down-regulated in S genotype; G4, 131 DEGs were up-regulated in S but down-regulated in R genotype; G5, 2,287 DEGs were only up-regulated in R genotype; G6, 1,843 DEGs were only down-regulated in R genotype; G7, 1,319 DEGs were only up-regulated in S genotype; G8, 1,781 DEGs were only down-regulated in S genotype.

    • The DEGs were assigned to enriched GO terms to annotate the gene function in terms of molecular function, cellular component, and biological process. Within cultivar expression differences, CM vs NM and CV vs NV, were the primary comparisons of interest. For up-regulated genes in CM vs NM, the major categories within molecular function were nucleoside phosphate binding (GO:1901265), carbohydrate derivative binding (GO:0097367) and purine nucleotide binding (GO:0017076). The top GO terms in cellular component were intracellular (GO:0005622) and membrane (GO:0016020), while single-organism process (GO:0044699), single-organism cellular process (GO:0044763) and biosynthetic process (GO:0009058) were identified as major biological process. For down-regulated genes in 'Meyer', the enriched biological processes were photosynthesis (GO:0015979) and alpha-amino acid biosynthetic process (GO:1901607), while enriched cellular components were membrane protein complex (GO:0098796) and thylakoid (GO:0009579). In 'Victoria', the up-regulated genes were mainly assigned to single-organism process (GO:0044699) and biosynthetic process (GO:0009058) in terms of biological process, intracellular (GO:0005622) and membrane (GO:0016020) in terms of cellular component, and ion binding (GO:0043167) in terms of molecular function. The down-regulated genes were enriched in biological process of transmembrane transport (GO:0098656, GO:1903825, GO:0003333), cellular component of photosystem I (GO:0009522, GO:0009538) and molecular function of transmembrane transporter activity (GO:0008514, GO:0005342, GO:0046943). The GO analysis of up-regulated DEGs is shown in Fig. 2, and more detail can be found in Supplemental File 2.

      Figure 2. 

      Gene Ontology (GO) classification of the up-regulated DEGs identified in (a) CM vs NM and (b) CV vs NV comparisons. DEGs were annotated in three categories: biological process, cellular component, and molecular function.

    • DEGs above were annotated using the KEGG database. Among the up-regulated genes in 'Meyer', 23 significant pathways were identified, with metabolic pathways, biosynthesis of secondary metabolites, and ribosome being the top three. Among down-regulated genes in 'Meyer', 13 pathways were identified, and top three were metabolic pathways, biosynthesis of secondary metabolites, and carbon metabolism. In the S genotype 'Victoria', up-regulated DEGs were assigned to 16 pathways, in which the top two pathways were the same as the R genotype, followed by plant hormone signal transduction. Meanwhile, the down-regulated DEGs were only assigned to two pathways: metabolic pathways and N-Glycan biosynthesis. The KEGG analysis of up-regulated DEGs is shown in Fig. 3, and more detail can be found in Supplemental File 3.

      Figure 3. 

      KEGG pathway classification of the up-regulated DEGs identified in (a) CM vs NM and (b) CV vs NV comparisons.

    • By searching against the reference genome and transcript, we identified a total of 674 genes from QTL regions previously identified as associated with freezing and winter injury traits, while 168 of them were overlapped with above DEGs list (Fig. 4; Table 2; Supplemental File 4). The DEGs were identified in all the 11 QTL regions except in QTL5 on chromosome 7. There were 74 genes identified in QTL1, which was the one with the highest number of genes, 16 genes in QTL2, six genes in QTL3, 19 genes in QTL4, nine genes in QTL6, 15 genes in QTL7, 12 genes in QTL8, 10 genes in QTL9, six genes in QTL10, and one gene in QTL11. We found that a lot of these genes have been reported to be involved in plant stress response, such as Late embryogenesis abundant protein, Heat stress transcription factor, and other transcription factors including NAC, WRKY, MYB, and ethylene-responsive transcription factors. There were also several signal transduction related genes like CBL-interacting protein kinase, auxin-responsive protein, and ABSCISIC ACID-INSENSITIVE 5-like protein. More details of these genes are shown in Supplemental File 4.

      Figure 4. 

      Heat map of partial candidate acclimation response genes that were identified both in the QTL regions previously reported in Brown et al.[15] and the differentially expressed genes reported here. Color bar shows the log2(Fold Change) value, red color indicates up-regulation, blue color indicates down-regulation, and white color indicates no significant change.

      Table 2.  Partial list of candidate acclimation response genes that were identified both in the QTL regions previously reported in Brown et al.[15] and differential expression genes reported in this study. Full gene list and QTL information are provided in Supplemental File 4.

      QTL
      region
      Gene_idlog2FoldChange 'Meyer'log2FoldChange 'Victoria'Blast gene accessionBlast gene name
      QTL1Zjn_sc00044.1.g05080.1.sm.mk−1.2035NSQ9C5Q2AI5L3_ARATH ABSCISIC ACID-INSENSITIVE 5-like protein 3
      QTL1Zjn_sc00044.1.g04640.1.sm.mkhc2.78124.4685Q7X996CIPK2_ORYSJ CBL-interacting protein kinase 2
      QTL1Zjn_sc00044.1.g04620.1.am.mkNS1.7587Q7XIW5CIPKT_ORYSJ CBL-interacting protein kinase 29
      QTL1Zjn_sc00044.1.g03340.1.am.mk3.9076NSQ7XHZ0HFB4B_ORYSJ Heat stress transcription factor B-4b
      QTL1Zjn_sc00044.1.g05110.1.sm.mkNS3.422A3AHG5LEA17_ORYSJ Late embryogenesis abundant protein 17
      QTL1Zjn_sc00044.1.g05130.1.sm.mkhc3.3773NSQ10MB4MYB2_ORYSJ Transcription factor MYB2
      QTL1Zjn_sc00044.1.g04960.1.sm.mk5.4299NSA0SPJ6NAMB2_TRITD NAC transcription factor NAM-B2
      QTL1Zjn_sc00044.1.g04580.1.am.mkhc−1.0768NSP42736RAP23_ARATH Ethylene-responsive transcription factor RAP2-3
      QTL1Zjn_sc00044.1.g04450.1.sm.mkhcNS1.8656Q9SW70SRP_VITRI Stress-related protein
      QTL2Zjn_sc00036.1.g01730.1.am.mk3.16174.869Q6VBA4HFC1A_ORYSJ Heat stress transcription factor C-1a
      QTL2Zjn_sc00036.1.g01650.1.am.mkhc3.3333NSQ93WV4WRK71_ARATH WRKY transcription factor 71
      QTL6Zjn_sc00007.1.g01680.1.am.mk4.04152.6034Q6K846IAA9_ORYSJ Auxin-responsive protein IAA9
      QTL8Zjn_sc00004.1.g04500.1.sm.mkhc−2.3169NSQ7XBH4MYB4_ORYSJ Transcription factor MYB4
      QTL8Zjn_sc00004.1.g04560.1.sm.mkhc1.3851−1.3166P14717PAL1_ORYSJ Phenylalanine ammonia-lyase
      QTL8Zjn_sc00004.1.g04570.1.sm.mkhc7.72815.9998A2X7F7PAL2_ORYSI Phenylalanine ammonia-lyase
      QTL10Zjn_sc00122.1.g00870.1.am.mk−1.0261NSO80341EF102_ARATH Ethylene-responsive transcription factor 5
      *log2FoldChange, positive value indicates up-regulation after cold acclimation, negative value indicates down-regulation after cold acclimation, NS indicates non-significant.
    • We previously identified 62 protein spots differentially accumulated in abundance under cold acclimation using 'Meyer' and 'Victoria'[16]. Five of them were repeatedly identified in the above DEGs, including Late embryogenesis abundant protein, Zinc finger protein, MAP kinase, Phospholipase D, and Photosystem II protein. There were 14 DEGs that belong to the Late embryogenesis abundant protein gene family, 16 to the MAP kinase gene family, 17 to the Zinc finger protein gene family, nine to the Phospholipase D gene family, and 12 to the Photosystem II protein gene family (Fig. 5; Table 3; Supplemental File 5). The gene Zjn_sc00044.1.g05110.1.sm.mk, which encodes a Late embryogenesis abundant protein (LEA17), was also found to be present in the QTL1 region listed above.

      Figure 5. 

      Heat map of Late embryogenesis abundant protein (LEA) genes that were identified as differentially expressed after cold acclimation in Z. japonica, which were also differentially expressed in protein abundance in a previous study[16]. Color bar shows the log2(Fold Change) value, red color indicates up-regulation, blue color indicates down-regulation, and white color indicates no significant change.

      Table 3.  List of Late embryogenesis abundant protein (LEA) genes that were identified as differentially expressed after cold acclimation in Z. japonica, which were also differentially expressed in protein abundance in a previous study[16].

      Gene_idlog2FoldChange
      'Meyer'
      log2FoldChange
      'Victoria'
      Blast gene
      accession
      Blast gene name
      Protein spot 64[16]6.7NSLate embryogenesis abundant protein 3 (LEA3)
      Zjn_sc00027.1.g00580.1.am.mk−1.34351.2295P46518LEA14_GOSHI Late embryogenesis abundant protein Lea14-A
      Zjn_sc00011.1.g01110.1.sm.mkNSInfQ94JF2LEA14_ORYSJ Late embryogenesis abundant protein 14
      Zjn_sc00011.1.g01130.1.sm.mk11.6469.6458Q94JF2LEA14_ORYSJ Late embryogenesis abundant protein 14
      Zjn_sc00044.1.g05110.1.sm.mkNS3.422A3AHG5LEA17_ORYSJ Late embryogenesis abundant protein 17
      Zjn_sc00002.1.g03150.1.am.mk3.80762.7007A3AHG5LEA17_ORYSJ Late embryogenesis abundant protein 17
      Zjn_sc00003.1.g02040.1.sm.mkhc7.5978.8459A3AHG5LEA17_ORYSJ Late embryogenesis abundant protein 17
      Zjn_sc00004.1.g08640.1.sm.mkInf9.7227Q96273LEA18_ARATH Late embryogenesis abundant protein 18
      Zjn_sc00027.1.g03080.1.am.mkhc4.28474.8089Q42376LEA3_MAIZE Late embryogenesis abundant protein, group 3
      Zjn_sc00012.1.g02810.1.sm.mkhcInf8.9382Q42376LEA3_MAIZE Late embryogenesis abundant protein, group 3
      Zjn_sc00081.1.g00290.1.sm.mk6.50693.7991Q9LJ97LEA31_ARATH Late embryogenesis abundant protein 31
      Zjn_sc00003.1.g11730.1.am.mkNSInfP09444LEA34_GOSHI Late embryogenesis abundant protein D-34
      Zjn_sc00108.1.g00280.1.sm.mkInf4.7263P09444LEA34_GOSHI Late embryogenesis abundant protein D-34
      Zjn_sc00023.1.g01980.1.am.mk−2.0172NSF4KFM8LEA65_ARATH Late embryogenesis abundant protein At5g17165
      Zjn_sc00004.1.g08970.1.am.mk−1.2806NSF4KFM8LEA65_ARATH Late embryogenesis abundant protein At5g17165
      * log2FoldChange, positive value indicates up-regulation after cold acclimation, negative value indicates down-regulation after cold acclimation, NS indicates non-significant, Inf indicates infinite value.
    • Cold acclimation, a low but above freezing temperature treatment, is a critical process to the development of freezing tolerance. The process directly inhibits metabolic reactions, and indirectly inhibits cold-induced osmotic stress, oxidative stress as well as numerous other genes and pathways by affecting lipid composition of membranes, sugars, and proline molecules, and increasing total soluble proteins[13,27]. To further disentangle the differences in metabolic response to freezing, the present study compared the transcriptomes of 'Meyer', a freeze-tolerant cultivar, and 'Victoria', a freeze-susceptible cultivar, under cold acclimated and non-cold acclimated conditions. The study identified abundant genes and pathways in response to cold acclimation, which provides insights into the molecular mechanisms underlying cold acclimation response in Z. japonica. In addition, the findings provide an opportunity to validate our previous studies on QTL mapping and proteomics, with the overlapped genes showing great value to be used for selection in breeding for freeze tolerance in Z. japonica.

    • As mentioned above, our previous studies have demonstrated that 'Meyer' and 'Victoria' perform differently when exposed to cold acclimation and freezing[1416]. The findings of this study provide evidence to explain these differences from the gene regulation mechanisms activated in response to cold acclimation. Cold response in plants is a very complex trait, which involves many different gene regulation and metabolic pathways[28]. The molecular basis of cold acclimation and acquired freezing tolerance in plants has been extensively studied, especially in model plant Arabidopsis. The CBF/DREB1-dependent cold signaling pathway has been characterized as the key cold signaling pathway[8,9]. In Arabidopsis, several CBF/DREB1s have been reported to be involved in the regulation of COR gene expression and cold tolerance[8,9]. In addition, the CBF/DREB1 (mainly CBF3/DREB1A) pathway was found to be controlled by transcription factor ICE1 (inducer of CBF expression1)[11]. In zoysiagrass, six DREB1 and two DREB2 unigenes were up-regulated significantly at the 2h-cold time point, and one DREB2 unigene was up-regulated at the 72h-cold time point using 'Meyer'[10]. DREB and ICE transcription factors were also identified in our study. There were two ICE1 transcription factors were found, Zjn_sc00034.1.g02330.1.am.mkhc was up-regulated in both cultivars, while Zjn_sc00094.1.g00080.1.am.mk was found in the G3 group (up-regulated in 'Meyer' but down-regulated in 'Victoria'), which indicates it might contribute to the higher cold tolerance in 'Meyer' (Supplemental File 1). There were 17 DREB transcription factors found, two of them in G1, two in G2, two in G5, six in G6, one in G7, and one in G8 (Supplemental File 1). Despite the common change trend in G1 and G2, more DREB genes were changed in 'Meyer' (G5, G6) than in 'Victoria' (G7, G8), indicating the two cultivars might have adopted different DREB members in response to cold stress.

      Expression of CORs (cold-regulated genes) regulated by CBF/DREB1s has been shown to be critical in plants for cold and freezing stress[29]. The cold-responsive protein kinase 1 (CRPK1), a plasma membrane protein which could phosphorylate 14-3-3 proteins, was proved to be a negative regulator of freezing tolerance. The crpk1 Arabidopsis mutant showed enhanced freezing tolerance[30]. There were two CRPK1 genes found in our study, Zjn_sc00002.1.g07100.1.sm.mkhc was up-regulated in both cultivars, while Zjn_sc00068.1.g00990.1.sm.mkhc was only down-regulated in 'Meyer' (Supplemental File 1). Considering CRPK1 gene is a negative regulator, decreased expression of Zjn_sc00068.1.g00990.1.sm.mkhc might enhance cold tolerance in 'Meyer'.

      To alleviate the oxidative damage in plant cells under cold stress, plants exhibit enhanced activities of antioxidant enzymes, such as superoxide dismutase (SOD), peroxidase (POD), catalase (CAT), glutathione peroxidase (GPX), and ascorbic acid peroxidase (APX) and eliminate excess reactive oxygen species (ROS)[31]. In this study, one, two, three and 34 genes were annotated as GPX, CAT, SOD and POD. The GPX gene was up-regulated in both cultivars. The two CAT genes were only changed in 'Victoria', one up-regulated and one down-regulated. The two SOD genes were up-regulated in both cultivars, while one SOD gene was only up-regulated in 'Meyer'. Five POD genes were down-regulated in both cultivars, three genes were up-regulated in both cultivars, three genes were up-regulated in 'Meyer' but down-regulated in 'Victoria', six genes only up-regulated in 'Meyer', four genes were only down-regulated in 'Meyer', one gene was only up-regulated in 'Victoria', while 12 POD genes were only down-regulated in 'Victoria' (Supplemental File 1). These findings indicated that more antioxidant genes were up-regulated in 'Meyer', and down-regulated in 'Victoria', suggesting that 'Victoria' might suffer more severe oxidative damage than 'Meyer' during cold stress.

    • Integrating QTL information and gene expression variation is a common strategy for candidate gene screening[3234]. In previous studies, we have identified multiple QTL associated with field winter injury, and lab-based post-freeze surviving green tissue and regrowth using a 'Meyer x Victoria' mapping population, and 11 QTL regions of interest were found that overlapped across different traits and under different environments. By searching genome sequences of these regions and integrating the DEGs identified in this study, we found a total of 168 genes spread across 10 of 11 QTL regions (Table 2; Supplemental File 4). Some of these have been proved to be involved in response to cold acclimation and freezing stress. Two CIPK (CBL-interacting protein kinase) genes were found, one was up-regulated in both cultivars, another one was only up-regulated in 'Victoria'. Arabidopsis CIPK3 and CIPK7 were reported to be involved in response to cold stress[35,36]. Zjn_sc00044.1.g04610.1.am.mk, annotated as a POD gene discussed above, was only down-regulated in 'Victoria', which might impair the freezing tolerance of 'Victoria'. Heat stress transcription factor (HSF) act as components of signal transduction that induce the expression of heat shock protein (HSP) in response to a broad range of abiotic stresses[37]. In this study, two HSF genes and a HSP gene located in QTL regions were found up-regulated by cold acclimation, while one of HSF gene, Zjn_sc00044.1.g03340.1.am.mk, only changed in 'Meyer' (Table 2). In addition, there were several transcription factors found in QTL regions, including MYB, NAC, ERF, WRKY, which function in regulating downstream cold response genes, and trigger a series of physiological and biochemical responses.

      Several genes were found overlapped with our previous comparative proteomic study between 'Meyer' and 'Victoria' under cold acclimation, including LEA, Zinc finger protein, MPK, Photosystem II (PSII) protein and Phospholipase (Table 3, Supplemental File 5). LEA proteins were identified in different plant species and proved that they are important factors in regulating plant freezing tolerance[12,38]. There were 14 LEA genes identified in DEGs, most of them were up-regulated in both cultivars, while three of them were down-regulated in 'Meyer'. It was noted that one LEA gene, Zjn_sc00044.1.g05110.1.sm.mk, was also found in the above QTL region, which shows great value to be used for cold tolerance improvement in Z. japonica. In plants, the cold response involves decreased photosynthetic capacity, and inhibition of photosystem II (PSII) activity is one of the most notable indicators of cold stress[39]. Recently, Zhang et al. reported that the expression of several PSII related genes were altered in Zoysia japonica during cold stress[40]. We identified a PSII protein in the proteomics study, which shows higher abundance in 'Meyer' but lower abundance in 'Victoria'[16]. In the present study, 12 genes encoding the PSII complex protein were found to be differentially expressed after cold acclimation. Most of the PSII genes were decreased during cold stress, while four of them were up-regulated in 'Victoria', and only one of them was up-regulated in 'Meyer'. The inconsistency between protein and transcript level might be due to post-translational modifications. Mitogen-activated protein kinases (MPK) were identified in both proteomics and transcriptomics levels in response to cold acclimation. MPK cascades are important signaling modules that convert environmental stimuli into cellular responses. Recent studies show that MPK3, MPK4, and MPK6 are rapidly activated after cold treatment[41]. Zhao et al. discovered that the MKK5-MPK3/MPK6 cascade negatively regulates cold response by promoting the degradation of ICE[42]. In contrast, the MEKK1-MKK2-MPK4 cascade positively regulates cold response, indicating the different roles of two MAP kinase cascades in the regulation of COR gene expression and freezing tolerance[43]. In this study, there were 16 MPK genes in DEGs, five of them were up-regulated and one was down-regulated in both cultivars. Expression of the rest of MPK genes showed different trends: three of them were only down-regulated in 'Meyer', and two genes were only down-regulated and five genes were only up-regulated in 'Victoria' (Supplemental File 5). These results indicated several different MPK pathways function in cold stress response in Z. japonica. Further studies are, however, necessary to demonstrate the function of these overlapped genes in Z. japonica response to cold acclimation and freezing stress.

    • This study examined the differentially expressed genes between freeze-susceptible 'Victoria' and freeze-tolerant 'Meyer' under cold acclimated and non-cold acclimated conditions. The differentially expressed genes showed different patterns in the two cultivars in response to cold acclimation. Integrating with previous studies, several genes were identified across the transcriptome, proteome, and QTL mapping to be associated with cold acclimation and freezing tolerance, such as LEA, CIPK, POD, HSF, HSP, MPK, PSII as well as multiple transcription factors. These findings strengthen the connection of these genes with resistance to cold and freeze stress, and upon further study, may be utilized to facilitate future breeding efforts to increase freeze tolerance in Z. japonica.

    • The authors confirm contribution to the paper as follows: funding acquisition: Milla-Lewis SR; study conception and design: Yu X, Milla-Lewis SR; data collection: Brown JM, Holloway HMP; resource and methodology: Tuong TD, Patton AJ, DaCosta M, Livingston DP; analysis and interpretation of results: Brown JM, Weldt CE, Yu X; draft manuscript preparation: Brown JM, Yu X; review and editing: Yu X, Milla-Lewis SR. All authors reviewed the results and approved the final version of the manuscript.

    • All the raw sequencing data were deposited into the database of sequence read archive (SRA) in National Center for Biotechnology Information (NCBI) under BioProject ID: PRJNA1040203. The other datasets generated during and/or analyzed during the current study are available from the corresponding author on reasonable request.

      • This research was supported in part by funding provided by the NC State University Plant Breeding Consortium, the NC State University Center for Turfgrass Environmental Research and Education, and the United States Golf Association. The authors are appreciative to Dr. Rui Shi from North Carolina State University for manuscript reviewing.

      • The authors declare that they have no conflict of interest. Susana R. Milla-Lewis is the Editorial Board member of Grass Research who was blinded from reviewing or making decisions on the manuscript. The article was subject to the journal's standard procedures, with peer-review handled independently of this Editorial Board member and the research groups.

      • Supplemental File 1 List of differentially expressed genes (DEGs) in different comparisons.
      • Supplemental File 2 Gene Ontology (GO) classification of the DEGs identified in different comparisons.
      • Supplemental File 3 KEGG classification of the DEGs identified in different comparisons.
      • Supplemental File 4 List of candidate acclimation response genes that be identified both in QTL regions and differentially expressed genes. QTL regions were previously reported in Brown et al.[15].
      • Supplemental File 5 List of candidate acclimation response genes that were also differentially expressed in protein abundance in a previous study[16].
      • Copyright: © 2023 by the author(s). Published by 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 (5)  Table (3) References (43)
  • About this article
    Cite this article
    Brown JM, Weldt CE, Holloway HMP, Tuong TD, Patton AJ, et al. 2023. Transcriptomic analysis of zoysiagrass (Zoysia japonica) provides novel insights into the molecular basis of cold acclimation. Grass Research 3:25 doi: 10.48130/GR-2023-0025
    Brown JM, Weldt CE, Holloway HMP, Tuong TD, Patton AJ, et al. 2023. Transcriptomic analysis of zoysiagrass (Zoysia japonica) provides novel insights into the molecular basis of cold acclimation. Grass Research 3:25 doi: 10.48130/GR-2023-0025

Catalog

  • About this article

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return