Phenotypic characterization and stem ultrastructure analysis of erect- and zigzag-type tea plants
Under natural conditions, the trees of MZ, QQ and LYQQ can grow upward uniformly. The leaves of MZ and QQ are flat, while those of LYQQ are folded inwards (Fig. S1). In the MZ plant, the stems grew straight up, exhibiting normal shoot morphology; however, the shoots of both QQ and LYQQ tended to bend at each node and were elongated in a zigzag fashion (Fig. 1A-C). To precisely investigate the differences in zigzag stems at the ultrastructure level, we dissected the longitudinal stems of the QQ, LYQQ and MZ tea plants. Observation of the stem longitudinal sections showed that the main tissue was basically normal, but the cell arrangement and shape differed between zigzag and erect stems (Fig. 1D-F). In QQ and LYQQ, the cortex cells tended to be disordered and arranged loosely, and the cells in both the cortex and pith exhibited an aberrant shape (Fig. 1G-I).
RNA sequencing, reference genome alignment and new gene annotation
To investigate the regulation of zigzag-stem formation in tea plants at the transcriptional level, we utilized RNA-Seq technology to analyse DEGs in the stems from MZ, QQ and LYQQ plants. In total, 69.12 Gb of raw data were generated from nine samples, and the sequence data were deposited in the NCBI Sequence Short Read Archive (SRA accession: PRJNA559220). After removing adaptor sequences, duplicate sequences, ambiguous reads and low-quality reads, a total of 24.16, 19.51, and 22.79 Gb clean reads were generated for MZ, QQ and LYQQ, respectively (Table S2). The average amount of clean reads per sample was 7.38 Gb. The Q20 values ranged from 97.37% to 98.55%, and the Q30 values ranged from 92.26% to 9507%. All of the transcripts were aligned to the reference genome, and the average proportion of the sample mapped to the genome was 76.38%. The new genes were then aligned to the Nr and KEGG databases for protein functional annotation. In total, 34248, 34374 and 33598 genes were identified from MZ, QQ and LYQQ, respectively, and 28021 (82.58%), 27441 (80.87%) and 27982 (82.46%) of genes were annotated as known genes in MZ, QQ and LYQQ, respectively (Table S2). These results indicated that the obtained high-quality transcriptomic data could be used for further analysis.
Validation of differential expression data
To validate the reliability of the RNA-Seq results, 16 DEGs were randomly selected from RNA-Seq data and examined using qRT-PCR. The qRT-PCR data exhibited similar expression patterns with the RNA-Seq data among cultivars (Fig. 2), suggesting that our transcriptomic data are reliable and valid.
Identification of DEGs and pathways between the comparisons of cultivars
The DEGs in each cultivar pair were then determined according to the parameters p value≤0.01 and |log2FC|≥1. In total, 6232 DEGs, including 2969 upregulated and 3263 downregulated DEGs, were detected in MZ-vs-QQ (Fig. 3A). GO enrichment analyses showed that these DEGs were classified into three categories: biological process, cellular component, and molecular function, and most of the DEGs were enriched in the terms ‘catalytic activity’, ‘metabolic process’, ‘cellular process’, ‘binding’, ‘single-organism process’, and ‘membrane’ (Fig. S2A). These DEGs were subjected to KEGG pathway enrichment analyses, which showed that the pathways ‘Biosynthesis of secondary metabolites’, ‘Plant-pathogen interaction’, ‘Phenylpropanoid biosynthesis’, ‘Stilbenoid, diarylheptanoid and gingerol biosynthesis’, ‘Monoterpenoid biosynthesis’, ‘Biosynthesis of unsaturated fatty acids’, ‘alpha-Linolenic acid metabolism’, and ‘Flavonoid biosynthesis’ were significantly enriched. Additionally, we found that the zeatin biosynthesis pathway was also enriched (Fig. 3B). In the MZ-vs-LYQQ comparison, a relatively high number of DEGs (7212), including 4002 upregulated and 3210 downregulated DEGs, were identified (Fig. 3A). Similar to the DEGs in MZ-vs-QQ, most of these DEGs were enriched in the terms ‘catalytic activity’, ‘metabolic process’, ‘cellular process’, ‘binding’, ‘single-organism process’, and ‘membrane’ (Fig. S2B). All the DEGs could be mapped to 132 KEGG pathways, and the pathways ‘Phenylpropanoid biosynthesis’, ‘Cutin, suberine and wax biosynthesis’, ‘Plant-pathogen interaction’, ‘Stilbenoid, diarylheptanoid and gingerol biosynthesis’, ‘Flavonoid biosynthesis’, ‘Biosynthesis of secondary metabolites’, ‘Monoterpenoid biosynthesis’, ‘Glutathione metabolism’, and ‘Arginine and proline metabolism’ were significantly enriched (Fig. 3B). A total of 6930 DEGs, including 3932 upregulated and 2998 downregulated DEGs, were detected in the QQ-vs-LYQQ comparison (Fig. 3A) and mapped to 132 pathways. The DEGs in the pathways ‘Plant-pathogen interaction’, ‘Cutin, suberine and wax biosynthesis’, ‘Biosynthesis of secondary metabolites’, ‘Phenylpropanoid biosynthesis’, ‘Stilbenoid, diarylheptanoid and gingerol biosynthesis’, ‘Brassinosteroid biosynthesis’ and ‘Monoterpenoid biosynthesis’ were significantly enriched (Fig. 3B).
Identification of DEGs and pathways involved in zigzag-stem formation in tea plants
We generated a Venn diagram to compare the different cultivars and showed that 1082 DEGs overlapped among the MZ-vs-LYQQ, MZ-vs-QQ, and QQ-vs-LYQQ comparisons (Fig. 4A). These DEGs were significantly enriched in the pathways of “Plant-pathogen interaction”, “Stilbenoid, diarylheptanoid and gingerol biosynthesis”, “Phenylalanine metabolism”, and “Tryptophan metabolism” (Fig. S3A). In addition, a total of 1255 DEGs, including 527 downregulated and 728 upregulated DEGs, were specifically detected in the MZ-vs-LYQQ comparison (Fig. 4A). Among the top 20 pathways, “Cysteine and methionine metabolism”, “Cutin, suberine and wax biosynthesis”, “Taurine and hypotaurine metabolism”, and “Other types of O-glycan biosynthesis” were markedly enriched (Fig. S3B). Unexpectedly, the number of DEGs in the MZ-vs-QQ set (949, including 494 downregulated and 455 upregulated) was lower than that in MZ-vs-LYQQ (Fig. 4A), and the pathways of “Glycosphingolipid biosynthesis − globo series” and “Limonene and pinene degradation” were significantly enriched (Fig. S3C). Additionally, a total of 1122 DEGs, including 593 upregulated and 529 downregulated DEGs, were specifically expressed in QQ-vs-LYQQ (Fig. 4A), but only the “Plant-pathogen interaction” pathway was considerably enriched (Fig. S3D). Moreover, a total of 2175 DEGs, including 1177 downregulated and 998 upregulated DEGs, overlapped between MZ-vs-LYQQ and MZ-vs-QQ specifically, indicating that these DEGs might be associated with zigzag-stem formation in both QQ and LYQQ. KEGG analysis showed that these DEGs were mainly involved in the “Plant-pathogen interaction”, “Phenylpropanoid biosynthesis”, “Flavonoid biosynthesis” and “Linoleic acid metabolism” pathways (Fig. 4B). GO enrichment analysis showed that these DEGs were significantly enriched in 59 GO terms, of which the most enriched components were categorized as catalytic activity (465), metabolic process (381), binding (323), cellular process (322), single-organism process (285) and membrane (274) (Fig. 4C).
Identification of key DEGs regulating zigzag-stem formation
Based on the changes in expression in the comparison of MZ-vs-QQ and MZ-vs-LYQQ, 76 DEGs potentially involved in zigzag-stem formation were identified (Fig. 5). Among these DEGs, 19 belonged to cell wall synthesis and cell expansion, of which seven, namely, cellulose synthase (TEA032164.1, TEA030545.1), expansin (TEA027164.1), leucine-rich repeat extensin-like protein 1 (XLOC_003301), chitinase-like protein (TEA022978.1) and pectinesterase (XLOC_003301, and TEA004581.1), were upregulated, whereas 12, especially xyloglucan endotransglucosylase/hydrolase (XLOC_007313, TEA019643.1, TEA031643.1), pectinesterase (TEA026842.1), reduced wall acetylation 2 (XLOC_021264), expansin (TEA012391.1), UDP-glycosyltransferase 74B5 (TEA020219.1) and isoamylase 3 (XLOC_040461), were downregulated in both QQ and LYQQ (Fig. 5A and Table S3). Eighteen transcription factor genes, including the 13 downregulated genes floral homeotic protein APETALA 1 (TEA017728.1), TIFY (TEA012041.1), NAC transcription factor 010 (TEA026168.1), transcription factor APETALA2 (XLOC_053049), WUSCHEL-related homeobox 2 (TEA032867.1), ethylene-responsive transcription factor TINY (TEA027175.1), transcription factor MYB1R1 (TEA026206.1), squamosa promoter-binding-like protein (TEA003577.1), transcription factor HEC1 (TEA030941.1), transcription factor bHLH18 (TEA000681.1), transcription factor SPATULA (TEA006216.1), growth-regulating factor 1 (TEA022970.1), and Scarecrow-like protein (TEA030046.1) and the five upregulated genes transcription factor bHLH041 (TEA031877.1), transcription factor IBH1-like (TEA009726.1), WRKY transcription factor 28 (TEA023233.1), transcription factor DIVARICATA (TEA031729.1) and transcription factor JUNGBRUNNEN 1 (TEA022287.1), were identified (Fig. 5B and Table S3). In addition, 10 DEGs involved in auxin, jasmonic acid, and salicylic acid metabolism and transport were also identified in the list of key DEGs; interestingly, except for jasmonic acid-amido synthetase (TEA020186.1), the remaining genes, especially PIN3 (TEA019069.1), were downregulated in both QQ and LYQQ (Fig. 5C and Table S3). Furthermore, seven DEGs involved in protein processing and transportation on the endoplasmic reticulum and vesicles, namely, vesicle-associated membrane protein 714 (XLOC_031693), SEC1 family transport protein, signal peptidase complex catalytic subunit SEC11A (TEA001395.1), SEC13B (XLOC_004426), SECA2 (XLOC_057225), SEC6 (TEA030236.1), SEC11A (TEA001395.1), and SEC22 (XLOC_037235); the three vacuolar protein sorting-related proteins VPS18 (TEA007337.1), VPS41 (TEA031089.1) and VSR6 (TEA021222.1); and the vacuole membrane protein KMS1 (XLOC_036914) were identified from the MZ-vs-QQ and MZ-vs-LYQQ sets (Fig. 5D and Table S3). Among these DEGs, VPS18, VPS41, SEC11A and SEC1 were significantly repressed in both QQ and LYQQ. The genes that regulate cell division (cell division cycle 20.1 and cell division control protein 6 B) and plant development, such as shaggy-related protein kinase, DEFECTIVE IN MERISTEM SILENCING 3 (XLOC_028596), RETICULATA-RELATED 3 (XLOC_032980), TOPLESS-like (XLOC_028345), TOPLESS-related protein (TEA008751.1), LAZY protein (TEA031847.1) and LAZY 1-like (TEA001744.1), were identified, and all of these genes were downregulated in both QQ and LYQQ (Fig. 5E and Table S3). Moreover, the DEGs VILLIN2 protein (VLN2) and actin-depolymerizing factor 2 (ADF2) were also suppressed in both QQ and LYQQ (Fig. 5E and Table S3).
Metabolic analysis and key metabolite identification
To investigate the metabolic pathways involved in zigzag-shaped stem formation, the metabolites in the stems of QQ, LYQQ and MZ were detected using UPLC-ESI-TOF-MS/MS. In total, 752 metabolites clustered into 97 KEGG pathways were identified from QQ, LYQQ and MZ, and among these metabolites, 75, 84 and 86 metabolites showed significantly different levels in the MZ-vs-QQ, MZ-vs-LYQQ and QQ-vs-LYQQ comparisons, respectively (Fig. S4). The Venn diagram analysis showed that 13 metabolites overlapped between MZ-vs-QQ and MZ-vs-LYQQ, which were our metabolites of interest (Fig. 6A), and the results indicated that these differential metabolites might be associated with zigzag-shaped stem formation in tea plants. Based on the log2 values of the fold changes, these differential metabolites were visualized as a heat map in Fig. 6B. Quercetin O-acetylhexoside, methyl gallate, D-pantothenic acid and L-glutamic acid were upregulated in both QQ and LYQQ, whereas the remaining metabolites, including fustin, 10-formyl-THF, skimmin, LysoPC 20:4, LysoPC 18:1 (2n isomer), LysoPC 18:3 (2n isomer), 2-methylsuccinic acid, 2-isopropylmalate, and caffeine, were significantly downregulated in tea plants with zigzag shoots.