-
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.
Sample Raw reads Clean reads Clean bases Error rate (%) Q20 (%) Q30 (%) GC content (%) NM1 49,475,254 48,077,382 7.2G 0.03 97.92 94.15 57.25 NM2 51,183,836 50,010,472 7.5G 0.03 97.84 93.99 57.49 NM3 50,168,628 49,079,396 7.4G 0.03 97.88 94.04 57.49 CM1 59,804,056 58,665,446 8.8G 0.03 97.37 92.95 52.47 CM2 67,204,258 65,599,380 9.8G 0.03 97.35 92.94 51.48 CM3 53,728,046 52,608,814 7.9G 0.03 97.11 92.48 52.52 NV1 71,039,902 69,474,296 10.4G 0.03 97.47 93.19 53.68 NV2 46,536,824 45,397,712 6.8G 0.03 97.39 93.03 52.99 NV3 55,073,384 53,918,022 8.1G 0.03 97.08 92.41 53.03 CV1 61,931,458 6,0855,104 9.1G 0.03 97.43 93.12 52.9 CV2 51,018,634 49,819,784 7.5G 0.03 97.28 92.84 52.43 CV3 67,843,324 66,500,252 10.0G 0.03 97.49 93.23 52.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. Differentially expressed genes (DEGs) analysis under acclimation treatment
-
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.
Gene Ontology (GO) classification of differentially expressed genes
-
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.
KEGG classification of differentially expressed genes
-
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.
Identification of acclimation response genes in the identified QTL regions
-
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
regionGene_id log2FoldChange 'Meyer' log2FoldChange 'Victoria' Blast gene accession Blast gene name QTL1 Zjn_sc00044.1.g05080.1.sm.mk −1.2035 NS Q9C5Q2 AI5L3_ARATH ABSCISIC ACID-INSENSITIVE 5-like protein 3 QTL1 Zjn_sc00044.1.g04640.1.sm.mkhc 2.7812 4.4685 Q7X996 CIPK2_ORYSJ CBL-interacting protein kinase 2 QTL1 Zjn_sc00044.1.g04620.1.am.mk NS 1.7587 Q7XIW5 CIPKT_ORYSJ CBL-interacting protein kinase 29 QTL1 Zjn_sc00044.1.g03340.1.am.mk 3.9076 NS Q7XHZ0 HFB4B_ORYSJ Heat stress transcription factor B-4b QTL1 Zjn_sc00044.1.g05110.1.sm.mk NS 3.422 A3AHG5 LEA17_ORYSJ Late embryogenesis abundant protein 17 QTL1 Zjn_sc00044.1.g05130.1.sm.mkhc 3.3773 NS Q10MB4 MYB2_ORYSJ Transcription factor MYB2 QTL1 Zjn_sc00044.1.g04960.1.sm.mk 5.4299 NS A0SPJ6 NAMB2_TRITD NAC transcription factor NAM-B2 QTL1 Zjn_sc00044.1.g04580.1.am.mkhc −1.0768 NS P42736 RAP23_ARATH Ethylene-responsive transcription factor RAP2-3 QTL1 Zjn_sc00044.1.g04450.1.sm.mkhc NS 1.8656 Q9SW70 SRP_VITRI Stress-related protein QTL2 Zjn_sc00036.1.g01730.1.am.mk 3.1617 4.869 Q6VBA4 HFC1A_ORYSJ Heat stress transcription factor C-1a QTL2 Zjn_sc00036.1.g01650.1.am.mkhc 3.3333 NS Q93WV4 WRK71_ARATH WRKY transcription factor 71 QTL6 Zjn_sc00007.1.g01680.1.am.mk 4.0415 2.6034 Q6K846 IAA9_ORYSJ Auxin-responsive protein IAA9 QTL8 Zjn_sc00004.1.g04500.1.sm.mkhc −2.3169 NS Q7XBH4 MYB4_ORYSJ Transcription factor MYB4 QTL8 Zjn_sc00004.1.g04560.1.sm.mkhc 1.3851 −1.3166 P14717 PAL1_ORYSJ Phenylalanine ammonia-lyase QTL8 Zjn_sc00004.1.g04570.1.sm.mkhc 7.7281 5.9998 A2X7F7 PAL2_ORYSI Phenylalanine ammonia-lyase QTL10 Zjn_sc00122.1.g00870.1.am.mk −1.0261 NS O80341 EF102_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. Identification of acclimation response genes that overlapped with previous proteomics study
-
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_id log2FoldChange
'Meyer'log2FoldChange
'Victoria'Blast gene
accessionBlast gene name Protein spot 64[16] 6.7 NS Late embryogenesis abundant protein 3 (LEA3) Zjn_sc00027.1.g00580.1.am.mk −1.3435 1.2295 P46518 LEA14_GOSHI Late embryogenesis abundant protein Lea14-A Zjn_sc00011.1.g01110.1.sm.mk NS Inf Q94JF2 LEA14_ORYSJ Late embryogenesis abundant protein 14 Zjn_sc00011.1.g01130.1.sm.mk 11.646 9.6458 Q94JF2 LEA14_ORYSJ Late embryogenesis abundant protein 14 Zjn_sc00044.1.g05110.1.sm.mk NS 3.422 A3AHG5 LEA17_ORYSJ Late embryogenesis abundant protein 17 Zjn_sc00002.1.g03150.1.am.mk 3.8076 2.7007 A3AHG5 LEA17_ORYSJ Late embryogenesis abundant protein 17 Zjn_sc00003.1.g02040.1.sm.mkhc 7.597 8.8459 A3AHG5 LEA17_ORYSJ Late embryogenesis abundant protein 17 Zjn_sc00004.1.g08640.1.sm.mk Inf 9.7227 Q96273 LEA18_ARATH Late embryogenesis abundant protein 18 Zjn_sc00027.1.g03080.1.am.mkhc 4.2847 4.8089 Q42376 LEA3_MAIZE Late embryogenesis abundant protein, group 3 Zjn_sc00012.1.g02810.1.sm.mkhc Inf 8.9382 Q42376 LEA3_MAIZE Late embryogenesis abundant protein, group 3 Zjn_sc00081.1.g00290.1.sm.mk 6.5069 3.7991 Q9LJ97 LEA31_ARATH Late embryogenesis abundant protein 31 Zjn_sc00003.1.g11730.1.am.mk NS Inf P09444 LEA34_GOSHI Late embryogenesis abundant protein D-34 Zjn_sc00108.1.g00280.1.sm.mk Inf 4.7263 P09444 LEA34_GOSHI Late embryogenesis abundant protein D-34 Zjn_sc00023.1.g01980.1.am.mk −2.0172 NS F4KFM8 LEA65_ARATH Late embryogenesis abundant protein At5g17165 Zjn_sc00004.1.g08970.1.am.mk −1.2806 NS F4KFM8 LEA65_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. -
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.
-
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
Transcriptomic analysis of zoysiagrass (Zoysia japonica) provides novel insights into the molecular basis of cold acclimation
- Received: 27 September 2023
- Accepted: 17 November 2023
- Published online: 05 December 2023
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.
-
Key words:
- Z. japonica /
- Cold acclimation /
- Freeze tolerance /
- Transcriptome /
- Differentially expressed genes