Transcriptome and metabolite analyses provide insights into zigzag-stem formation in tea plants ( Camellia sinensis )

Shoot orientation is important for plant architecture formation, and the zigzag shoot is a special trait of many plants. The zigzag-shaped shoot has been selected and well studied in Arabidopsis , however, the regulatory mechanism of zigzag-shaped shoot development in other plants especially woody plants is largely unknown. In this study, tea plants with zigzag shoots, namely, Qiqu (QQ) and Lianyuanqiqu (LYQQ), were investigated and compared with the erect-shoot tea plant Meizhan (MZ) in an attempt to reveal the regulation of zigzag-shoot formation. Tissue section observation showed that the cell arrangement and shape of zigzag stems were aberrant compared with those of normal shoots. Moreover, a total of 2175 DEGs were identified from the zigzag shoots of the tea plants QQ and LYQQ compared to the shoots of MZ using transcriptome sequencing, and the DEGs involved in the “Plant-pathogen interaction”, “Phenylpropanoid biosynthesis”, “Flavonoid biosynthesis” and “Linoleic acid metabolism” pathways were significantly enriched. Additionally, the DEGs associated with cell expansion, vesicular trafficking, phytohormones, and transcription factors were identified and analysed. Metabolomic analysis showed that 13 metabolites overlapped and were significantly changed in the shoots of QQ and LYQQ compared to MZ. Our results suggest that zigzag-shaped shoot formation might be associated with the gravitropism response and polar auxin transport in tea plants and provides a valuable foundation for further understanding the regulation of plant architecture formation and for the cultivation and application of horticultural plants in the future. of the regulatory molecular mechanism of zigzag-shoot formation in woody plants. In the current study, we investigated the mechanism of formation in tea plants using comparative transcriptomics and metabolomic analysis. The results showed that zigzag-shaped stem development in tea plants might be regulated by a complex network involving vesicular trafficking, phytohormones, cell expansion, secondary metabolism, and diverse transcription factors. Importantly, zigzag-shoot formation might be closely related to alterations in the gravitropic response and polar auxin transport in tea plants. Our findings provide insights into zigzag-shaped shoot formation in tea plants and also serves as a valuable foundation for further investigations of architecture formation in woody plants.

Gene abundances were quantified by RSEM software, and the unigene expression levels were quantified using FPKM (fragments per kilobase of transcript per million mapped reads) values [30]. To identify DEGs across samples or groups, the edgeR package (http://www.rproject.org/) was used, and genes with |log2FC| ≥1 and FDR <0.05 were identified as DEGs. DEGs were then subjected to enrichment analysis of GO functions and KEGG pathways using the GOseq R package [31] and KOBAS software [32], respectively.

Quantitative real-time PCR validation
To validate the reliability of the gene expression analysis in this study, 15 genes were selected for real-time qPCR analysis. Total RNA samples were used for cDNA synthesis according to the method recommended by the kit manufacturer. The primer information for qRT-PCR analysis is listed in Table   S1. qRT-PCR was performed using SYBR Premix Ex Taq™ II (TaKaRa, Dalian, China) in a CFX96 Touch real-time PCR system (BIO-RAD, California, USA) according to the manufacturer's protocol, and the amplification was performed as our previously reported [33]. The results were calculated using the 2 −ΔΔCT method [34] with the CsPTB1 gene as a control [35]. Each sample was examined in triplicate.

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

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. To understand the cause of the zigzag morphological anomaly, we observed stem sections using optical microscopy. There were no aberrations in the main tissue structure (Fig. 1). However, differences in cell arrangement and shape were observed between the zigzag stems of QQ and LYQQ and the MZ stems, including the arrangement of cortex cells, which were disordered and loosely arranged, and the cells in the cortex and pith exhibited aberrant shapes (Fig. 1). Consistent with this finding, a similar cell shape and arrangement were also found in the zig-1 stems in Arabidopsis [8], indicating that the zigzag shape of the stems might be caused by the anomalous shape, arrangement and expansion of cells.
To investigate the molecular mechanism of zigzag-shaped stem formation in tea plants, the shoots were collected, and transcriptome sequence analysis was performed. In total, 69.12 Gb of sequences were generated, and 6232, 7212 and 6930 DEGs were identified from MZ-vs-QQ, MZ-vs-LYQQ and QQvs-LYQQ, respectively (Table S2 and Fig. 3B). These DEGs were mainly enriched in several pathways, such as 'Biosynthesis of secondary metabolites', 'Plant-pathogen interaction', and 'Stilbenoid, diarylheptanoid and gingerol biosynthesis', indicating that these pathways might be associated with differences among cultivars (Fig. 3B). To gain insights into the DEGs that specifically regulate zigzag stems, we made a Venn diagram for the cultivar comparisons and identified 2175 overlapping DEGs between MZ-vs-QQ and MZ-vs-LYQQ, which were mainly enriched in the "Plant-pathogen interaction", "Phenylpropanoid biosynthesis", "Flavonoid biosynthesis" and "Linoleic acid metabolism" pathways ( Fig. 4). Among these pathways, the "Phenylpropanoid biosynthesis" pathway serves as a source of degree inclined radiate pine seedlings was higher compared to the lower half and to non-inclined seedlings when seedlings were inclined for one month, indicated that quercetin can accumulate in the inclined shoots such as internodes of the zigzag shoots. Consequently, we observed that a C4H gene (TEA016772.1) was significantly expressed between MZ-vs-QQ and MZ-vs-LYQQ; moreover, we found that skimmin (pmf0295) expression was reduced in the MZ-vs-QQ (log2 FC: -0.995), and MZ-vs-LYQQ (log2 FC: -3.52) comparisons. Therefore, the zigzag-shaped stems of tea plants might be partially related to flavonoids, especially flavonols mediated auxin transport.
It is well recognized that the zigzag inflorescence stem is a result of the mutation of zig (zigzag)/sgr4, which encodes VPS10 interacting 11 (VTI11), a Qb-SNARE involved in vesicle transport between the trans-Golgi network and vacuole that exhibits abnormal gravitropism in Arabidopsis  Table S3). Among these genes, both vacuolar protein sortingassociated protein genes (VPS18 and VPS41) were repressed on QQ and LYQQ, whereas the expression levels of vesicle-associated membrane protein 714 (VAMP714) and vacuolar-sorting receptor 6 (VSR6) were upregulated in QQ and LYQQ. It has been established that the phenotype of zig-1 could be partially suppressed by mutation of the zig suppressor1 (ZIP1), ZIP2, ZIP3 and ZIP4 genes [46, 51, 52]. Niihama, et al. [52] reported that ZIP2, encoding an AtVPS41/AtVAM2 protein, is involved in protein sorting to vacuoles in Arabidopsis, and the zip2 mutation is a missense mutation.
These results indicate that zigzag-shaped stem formation is mainly related to abnormal gravitropism responses mediated by membrane trafficking. Additionally, we also found that six SEC family genes were significantly expressed in the MZ-vs-QQ and MZ-vs-LYQQ sets (Fig. 5D and Table S3), which is crucial for SNARE complex assembly and preprotein translocation [53,54]. Interestingly, we found that both LAZY genes (TEA031847.1, TEA001744.1) were markedly suppressed in QQ and LYQQ ( Fig.   5D and Table S3). It is well known that loss-of-function mutation of LAZY1 enhances polar indole-3acetic acid (IAA) transport and reduces shoot gravitropism and therefore regulates the growth angle of lateral branches in rice, Arabidopsis and maize [15-17, 19, 55]; however, the function of LAZY in woody plants remains to be elucidated. Therefore, we propose that the development of zigzag-shaped stems of tea plants might be associated with a change in the shoot gravitropism response, especially that affecting the disruption of membrane trafficking to vacuoles, although the zig/sgr4 gene did not exhibit a significant change.
The plant hormone auxin is important for organ growth and cell morphogenesis. In this study, seven DEGs involved in auxin metabolism, transport and signalling were identified by comparison of the MZvs-QQ and MZ-vs-LYQQ sets, and all of the genes were downregulated in QQ and LYQQ ( Fig. 5C and Table S3). Numerous studies have suggested that polar auxin transport, which is primarily determined  Table S3); these changes in expression in QQ and LYQQ might alter the polar transport of auxin and then affect the auxin gradients between stem sides; therefore, the shoot exhibits bending. Wu, et al.
[59] reported that VLN2, a type of actin-binding protein involved in microfilament regulation, affects the recycling of PIN2 and polar auxin transport, and vln2-defective mutant plants exhibited a hypersensitive gravitropic response and twisted roots and shoots at the seedling stage. In this study, we also found that VLN2 (TEA030753.1) was markedly repressed in QQ and LYQQ ( Fig. 5E and Table   S3), suggesting that the zigzag-shaped stems in tea plants might be related to polar auxin transport and the gravitropism response.
Moreover, cell expansion might exert a compressive force, leading to bending of the stem. Our results showed that the genes involved in cell expansion and cell wall synthesis, such as expansin, REDUCED WALL ACETYLATION, and xyloglucan endotransglucosylase/hydrolase protein, were differentially expressed ( Fig. 5A and Table S3), resulting in alteration of cell expansion and elongation. In addition, WAK genes were also differentially expressed in MZ-vs-QQ and MZ-vs-LYQQ, suggesting that WAK-mediated cell expansion and signalling pathways might be required for zigzag-shaped stem formation in tea plants. It is possible that cell expansion in the stems can produce a force that can lead to zigzag-shaped shoots. Importantly, cell expansion and differentiation predominantly rely on auxin [11, 60, 61]; thus, the mechanism of auxin regulation in zigzag-shoot development needs to be studied precisely in the different tissue sides of zigzag stems.
In this study, 20 TFs (including two ARFs) belonging to different TF families were differentially expressed in MZ-vs-QQ and MZ-vs-LYQQ, and most of these TFs were significantly repressed in QQ and LYQQ (Fig. 5B and Table S3). Almost all of the homologous genes were associated with plant Angle Control 1 (TAC1) and WEEP in peach [21, 67, 68], have been isolated; however, most of these genes are related to plant gravity response, and none of these genes encode TFs. Additionally, among the nine SGR genes, SGR1 and SGR7 encode scarecrow (SCR) and GRAS family TFs, respectively. We found that the SCR gene (TEA030046.1) was markedly repressed in QQ and LYQQ ( Fig. 5B and Table   S3). Plant growth and development is conducted by a TF-mediated complex network integrating with plant hormones, enzymes and other cellular component; therefore, as-yet-unknown TFs are possibly involved in the regulation of plant architectural phenotypes.

Conclusion
In the current study, we investigated the mechanism of zigzag-shoot formation in tea plants using

Supplementary Files
This is a list of supplementary files associated with this preprint. Click to download.