Phenotypic characterization and stem ultrastructure analysis of tea plants with erect and zigzag-shaped shoots
Under natural conditions, the trees of MZ, QQ and LYQQ can grow upward uniformly (Additional file 1: Fig. S1). The leaves of MZ and QQ are flat, while those of LYQQ are folded inwards. 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 elongate in a zigzag fashion (Fig. 1a-c). Additionally, plants with zigzag-shaped shoots had shorter stems and fewer leaves than erect plants (Additional file 1: Fig. S2a), and the internode distance (between two nodes) in plants with zigzag-shaped shoots was significantly shorter than that in erect plants (Additional file 1: Fig. S2b). To precisely investigate the differences in zigzag-shaped stems at the ultrastructure level, we longitudinally dissected the stems of the QQ, LYQQ and MZ tea plants. Observation of the stem longitudinal sections showed that the tissues were basically normal, but the cell arrangement and shape differed between zigzag-shaped 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 aberrant shapes (Fig. 1g-i).
RNA sequencing, reference genome alignment and new gene annotation
To investigate the regulation of zigzag-shaped stem formation in tea plants at the transcriptional level, we utilized RNA-Seq technology to analyse DEGs in the stems of MZ, QQ and LYQQ plants. In total, 46.06 million clean reads 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 16.37, 13.28, and 15.48 million high-quality clean reads were generated for MZ, QQ and LYQQ, respectively (Additional file 2: Table S1). The average amount of clean reads per sample was 5.2 million. The Q20 values ranged from 97.39% to 98.55%, and the Q30 values ranged from 92.26% to 95.07%. All the transcripts were aligned to the reference genome, and the average proportion of samples 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%) genes were annotated as known genes in MZ, QQ and LYQQ, respectively (Additional file 2: 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 the RNA-Seq data and examined using qRT-PCR. The qRT-PCR data exhibited similar expression patterns to the RNA-Seq data among the cultivars (Fig. 2), suggesting that our transcriptomic data are reliable and valid.
Identification of DEGs and pathways in cultivar comparisons
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 most of the DEGs were enriched in the terms ‘catalytic activity’, ‘metabolic process’, ‘cellular process’, ‘binding’, ‘single-organism process’, and ‘membrane’ (Additional file 1: Fig. S3). The DEGs were also 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 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). 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-shaped 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” (Additional file 1: Fig. S4a). 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 (Additional file 1: Fig. S4b). 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 “Glycosphingolipid biosynthesis − globo series” and “Limonene and pinene degradation” were significantly enriched (Additional file 1: Fig. S4c). 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 significantly enriched (Additional file 1: Fig. S4d). 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-shaped 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 highly 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-shaped stem formation
Based on the changes in expression in the comparisons MZ-vs-QQ and MZ-vs-LYQQ, 76 DEGs potentially involved in zigzag-shaped stem formation were identified (Fig. 5). Among these DEGs, 19 were associated with 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 Additional file 3: 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 Additional file 3: 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 Additional file 3: 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 comparisons (Fig. 5d and Additional file 3: Table S3). Among these DEGs, VPS18, VPS41, SEC11A and SEC1 were significantly repressed in both QQ and LYQQ. 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 also identified, and all of these genes were downregulated in both QQ and LYQQ (Fig. 5e and Additional file 3: Table S3). Moreover, the DEGs VILLIN2 protein (VLN2) and actin-depolymerizing factor 2 (ADF2) were suppressed in both QQ and LYQQ (Fig. 5e and Additional file 3: 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 (Additional file 1: Fig. S5). 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 their log2 fold change values, 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-shaped shoots.