Three tea varieties, Feiyun (FY), Qiancha 1 (QC1) and Longjing 43 (LJ43) were planted in tea germplasms in the Germplasm Tea Repository of the Guizhou Tea Research Institute located in Meitan County (N27°45’, E107°29’), Zunyi City, Guizhou Province, China. The reason for choosing these three varieties is that significant differences in stem length were found in the investigation of their stem length. In the spring of 2024, when the tea plants reached the stage of one bud and two leaves, stem tissues of the three tea plant varieties were sampled. The stem tissue between the second and fish leaves of these varieties was used as internode samples ( Figure 1A ), and three biological replicates were collected for each group. All samples were frozen in liquid nitrogen and stored in dry ice. They were brought back to the laboratory and stored at -80°C until use.

When FY, QC1, and LJ43 reach the stage of one bud and two leaves, take 15 branches each from bud to fish leaf as the samples. The representative samples of the three varieties were photographed and the differences in internode length were identified according to their phenotypes. In addition, using a ruler to accurately measure the internode length between the second leaf and the fish leaf of all samples. Record and statistically analyze the measurement results of all samples from three groups, with accuracy rounded to one decimal place.

Nine stem tissues from three varieties were ground in liquid nitrogen and total RNA was extracted using RNAprep Pure Plant Plus Kit (cat DP441, Tiangen, China), according to the manufacturer’s protocol. The quality and quantity of each RNA extract was determined using automated electrophoresis system (Fragment Analyzer 5400) and Qubit fluorescence quantitative analyzer. The qualified RNA was reverse-transcribed and the chain specific library was built for transcriptome sequencing in Lianchuan Biotechnology Co., LTD. Paired-end 150-bp sequencing raw reads were obtained. Clean reads were generated from next-generation sequencing data after removing adapter and primer sequences and low-quality reads using Cutadapt (Martin, 2011). High-quality clean reads were used for subsequent analyses. The raw RNA-seq data are publicly available at the Genome Sequence Archive database of China National Center for Bioinformation (https://ngdc.cncb.ac.cn/gsa) under project accession number PRJCA036838.

Splice-site and exon information were extracted from GTF file with Qiancha 1 genome (Qiao et al., 2025) annotation using Python scripts in HISAT2 (v.2.2.1) software. Then constructing a reference genome sequence index using the -build parameter of HISAT2 software. All clean data were mapped to the tea plant reference genome using HISAT2 with default parameter to generate SAM files. Samtools software was used to convert Sam files into BAM files (Ramirez-Gonzalez et al., 2012). Transcripts were assembled using the BAM files by StringTie software (v.2.2.1) with default settings (Pertea et al., 2016). To quantify the transcript expression level, the fragments per kilobase of per million mapped fragments (FPKM) were calculated by StringTie. To ensure the reliability of the assembly results, if a transcript has FPKM values less than one in three biological replicates of one sample, it will be removed and the remaining transcript will be considered to exist in the sample. The FPKM of the expressed genes were extracted in each sample and formed a table. The Non-redundant (Nr), Kyoto Encyclopedia of Genes and Genomes (KEGG), and Gene Ontology (GO) annotation information of the reference genome were downloaded for further analysis of the functions of these genes.

DEGseq package from R software (v4.2.2) was used to calculate differentially expressed genes (DEGs) based on their FPKM values in different groups (Wang et al., 2010). If the fold-change > 2 and q-values < 0.05 of a gene in two sets of samples, we define this gene as DEG. All DEGs between the two groups were statistically analyzed after removing redundancy and compared using a Venn diagram. TBtools software (Chen et al., 2023) was used for GO and KEGG analysis by their background files and DEGs. GO terms with corrected P < 0.01 were considered significantly enriched and visualizing the enrichment result by t R software.

The stem internode length and expression of DEGs were used for co-expression analysis by the WGCNA package of R software (Pei et al., 2017). Cluster analysis is performed on DEGs, and the weighted values of correlation coefficients are calculated to make the genes in the network follow a scale-free network distribution. Selecting appropriate power values and constructing a one-step network to divide differentially expressed genes into different modules. The detailed parameters are as follow: power = 24, TOMType = unsigned, and branch merge cut height = 0.25. Finally, the constructed gene network and internode length were correlated to identify the gene modules significantly related to the trait.

Total RNA was synthesized into the first strand cDNA using the PrimeScript RT Reagent Kit (cat RR036A, Takara, Japan) using the manufacturer’s protocols. The cDNA of all samples was adjusted to a consistent concentration. The 10ul PCR amplification includes 5 μ L of TB Green (cat RR820A, Takara, Japan), 0.2 μ L of upstream and downstream primers, 1.5 μ L cDNA template, and 3.1 μ L of sterile water, and thermal cycling is performed using a two-step method. All reactions were run in technical triplicates and the CsActin gene was selected as the internal control (Mi et al., 2025). Relative expression of selected genes were analyzed using the 2^-△Ct method. All primers used in qRT-PCR were listed in Supplementary Table S1 .

In order to observe the differences in stem internode length among these three tea plant varieties, we measured the internode length from the second leaf to the fish leaf during their one bud and two leaf stages ( Figure 1A ). The results showed that FY had the longest internode length, followed by QC1 and LJ43, with LJ43 having the shortest internode length. The statistical results of measuring 15 internode lengths of each of the three varieties were consistent with the photographic results. The internode length of stem in FY, QC1 and LJ43 showed significant difference at p < 0.05 ( Figure 1B ). FY internode length ranges from 1.8 to 3.8 cm, with an average length of 2.4 cm; QC1 internode length ranges from 1.4 to 2.6 cm, with an average length of 1.9 cm; LJ43 has the shortest internode length, ranging from 1.1 to 2.3 cm with an average length of 1.6 cm ( Figure 1B ; Supplementary Table S2 ). The internode length of FY is significantly greater than that of QC1 and LJ43, and there is also a difference between QC1 and LJ43, but less than the difference between FY and them.

To identify potential genes regulating stem length, we performed RNA-seq on stem tissues of FY, QC1, and LJ43. RNA-seq obtained a total of 48.85 Gb of raw data (range 5.32 ~ 6.31 Gb), Q30 was higher than 96%, GC content was about 44% ( Supplementary Table S3 ). The DEGs in the three groups were analyzed, and a total of 10,518 DEGs were identified. Among them, there are 2566 DEGs including 1602 upregulated and 964 downregulated genes in QC1 vs FY, 6703 DEGs including 3907 upregulated and 2796 downregulated genes in QC1 vs LJ43, and 7859 DEGs including 3336 upregulated and 4523 downregulated genes in FY vs LJ43 ( Figure 2A ). We found that the two varieties with the largest internode length difference (FY and LJ43) also had the largest number of DEGs, while QC1H vs FY had the smallest number of DEGs. Furthermore, we performed a Venn diagram analysis of the DEGs in the three comparison groups. The result shows that only 542 of the 10,518 DEGs were differentially expressed in all three comparison groups. 1135, 272, and 2225 DEGs were identified individually only in groups QC1H vs LJ43, QC1 vs FY, and FY vs LJ43 ( Figure 2B ).

Heat maps were used to analyze the expression patterns of DEGs in the three varieties. According to their expression patterns, they can be classified into three categories ( Figure 2C ). The expression level of the first type of DEGs are highest in FY and relatively low in QC1 and LJ43; The expression level of the second type of DEGs are highest in QC1 and relatively low in FY and LJ43; The DEGs in the third category are highest in LJ43, while they are relatively low in FY and QC1. In addition, some DEGs in the first category showed a gradually decreasing trend in FY, QC1, and LJ43, while some DEGs in the third category showed a gradually increasing trend in FY, QC1, and LJ43. The expression trends of these two differentially expressed genes show the same or opposite relationship with the internode length of their stems, and may play a positive or negative regulatory role in the internode length.

To further investigate the functions of these DEGs, we conducted GO and KEGG enrichment analyses on them ( Figure 3 ). The GO and KEGG enrichment pathways with p-values less than 0.01 after correction are displayed through bar charts and bubble plots. There are 5, 4, and 11 pathways for DEGs genes in the molecular function, cellular components, and biological processes modules of GO analysis. In molecular function, DEGs are mainly enriched in DNA-binding transcription factor activity, transcription regulator activity and DNA binding. In biological processes, DEGs are mainly enriched in signal transduction and biosynthetic process. KEGG enrichment analysis found that DEGs were mainly enriched in plant hormone signal transduction, signal transduction and flavonoid biosynthesis. It is worth noting that both KEGG and GO analyses enriched signal transduction pathways, indicating that hormone signals or signal receptor proteins may play an important role in regulating internode length. In addition, the regulation of transcription factors is also a major pathway for enrichment.

To explore the DEGs related to internode length of tea plants, we conducted WGCNA analysis on the DEGs identified in three tea varieties and their internode length. The results showed that 10,518 DEGs were clustered into 9 modules based on their expression patterns ( Figure 4A ), namely Black (106 DEGs), Blue (3168 DEGs), Pink (64 DEGs), Green (4412 DEGs), Yellow (441 DEGs), Red (253 DEGs), Brown (1679 DEGs), Turquoise (356 DEGs) and Grey (39 DEGs). In order to further identify DEGs related to stem internode length, we set a screening criterion and considered modules with | r | > 0.8 and P < 0.01 in the correlation analysis results as modules significantly correlated with internode length phenotype. According to this standard, we identified a blue module (r = 0.89, P = 0.001) significantly positively correlated with stem internode length and a green module (r = -0.95, P= 1e-4) significantly negatively correlated with this phenotype. Interestingly, DEGs positively and negatively correlated with internode length were identified in the 9 modules, indicating that gene regulation of tea stem internode length is influenced by both positive and negative factors.

Through the enrichment analysis of DEGs, it was found that signal transduction and transcription factors may play a role in regulating stem elongation. In addition, gene modules significantly positively and negatively correlated with internode length were identified through WGCNA analysis. To further screen for key candidate genes, we performed Venn diagram analysis on the WGCNA results and enrichment results. There are a total of 3168 DEGs in the blue module positively correlated with stem length, and 581 DEGs in the enriched transcription factor pathway, of which 246 are significantly positively correlated with stem length ( Figure 4B ). There are 256 DEGs enriched in the signal transduction pathway, and 154 DEGs intersect with the blue module. Using the same method, the DEGs negatively correlated with stem length in the green module were intersected with DEGs in transcription factors and signal transduction pathways, and 331 transcription factors negatively correlated with stem length and 181 signal transduction genes negatively correlated with stem length were identified. Using word cloud analysis to analyze these transcription factors related to stem length, it was found that these transcription factors are mainly ethylene−responsive transcription factor (ERR), WRKY transcription factor, transcription factor bHLH, transcription factor MYB, and NAC domain−containing protein ( Figure 4C ).

These identified candidate genes were used as nodes to construct a co expression network of genes related to stem internode length. In the blue module, a total of 22 candidate genes related to stem length were identified, and these candidate genes were classified into 6 categories according to their gene annotations. Among them, there were 5 types of transcription factors, namely 5 WRKY, 6 ERF, 3 NAC, 2 bHLH, and 3 MYB. In addition, 3 gibberellin 2-oxidase genes (GA2ox) in gibberellin metabolic pathways were identified ( Figure 5 ). Based on the number of candidate gene nodes, they are represented by circles of different sizes. The information of all candidate genes is listed in Supplementary Table S4 . For example, in the blue module, ERF (CsQC014983) has a maximum of 2683 connections with other genes, while MYB (CsQC027741) has a minimum of 13 connections with other genes ( Supplementary Table S4 ). In the green module, a total of 6 candidate genes related to stem length were identified, including 3 bHLH transcription factors, 1 WRKY transcription factor, 1 DELLA protein, and 1 gibberellin 3-oxidase genes (GA3ox) gene. A total of 28 hub genes were identified in two modules, which may play a positive or negative role in regulating tea stem elongation.

The expression patterns of candidate genes related to internode length in three tea varieties were analyzed by qRT-PCR. The results showed that the selected 9 genes were differentially expressed in all three varieties ( Figure 6 ). The expression level of one DELLA protein (CsQC006443) gradually increased in FY, QC1, and LJ43, and there were significant differences among the three varieties. This change pattern is negatively correlated with the internode length of the three varieties. Another DELLA protein (CsQC023224) is highly expressed in FY and QC1, while its expression level is lower in LJ43. The expression trends of two GA2ox genes (CsQC027290 and CsQC041233) in the three varieties are consistent, with high expression in FY, followed by QC1 with the lowest expression level in LJ43. The expression trend of GA3ox (CsQC053710) is opposite to that of GA2ox, with gradually increasing expression levels in FY, QC1, and LJ43. The expression patterns of MYB (CsQC027741), ERF (CsQC030079), and two WRKY transcription factors (CsQC050228 and CsQC053961) in the three varieties are consistent, with the highest expression in FY followed by QC1, and the lowest expression level in LJ43. The expression patterns of these four transcription factors in the three varieties are significantly positively correlated with their stem internode length, which may play a role in regulating stem elongation.