Comparative Microscopic, Transcriptome and IAA Content Analyses Reveal the Stem Growth Variations in Two Cultivars Ilex verticillata

Ilex verticillata is not only an excellent ornamental tree species for courtyards, but it is also a popular bonsai tree. ‘Oosterwijk’ and ‘Red sprite’ are two varieties of Ilex verticillata. The former has a long stem with few branches, while the latter has a short stem. In order to explain the stem growth differences between the two cultivars ‘Oosterwijk’ and ‘Red sprite’, determination of the microstructure, transcriptome sequence and IAA content was carried out. The results showed that the xylem thickness, vessel area and vessel number of ‘Oosterwijk’ were larger than in ‘Red sprite’. In addition, our analysis revealed that the differentially expressed genes which were enriched in phenylpropanoid biosynthesis; phenylalanine metabolism and phenylalanine, tyrosine and tryptophan biosynthesis in the black and tan modules of the two varieties. We found that AST, HCT and bHLH 94 may be key genes in the formation of shoot difference. Moreover, we found that the IAA content and auxin-related DEGs GH3.6, GH3, ATRP5, IAA27, SAUR36-like, GH3.6-like and AIP 10A5-like may play important roles in the formation of shoot differences. In summary, these results indicated that stem growth variations of ‘Oosterwijk’ and ‘Red sprite’ were associated with DEGs related to phenylpropanoid biosynthesis, phenylalanine metabolism and phenylalanine, tyrosine and tryptophan biosynthesis, as well as auxin content and DEGs related to the auxin signaling pathway.


Introduction
The whorled holly (Ilex verticillata), also known as north American holly, is an excellent ornamental fruit tree species. 'Oosterwijk', the most popular whorled holly variety in the cut branch market, has long branches which provide a good basis for the shape of the cut branches. Another variety, 'Red sprite', has a short plant type with many branches and is mainly used in the bonsai market. Plant height is a key agronomic trait closely related to plant structure, photosynthetic efficiency and productivity [1]. Height is determined by several developmental factors such as the number of phytomers and whether a plant genotypes, microstructural observations revealed different lignin contents between the two varieties and the stem lignin synthesis gene was screened by transcriptome data [28]. In Cynodon dactylon L., microstructure and transcriptome sequencing were used to study the genetic control of the tiller angle and stem growth [29]. The differences in stem bending between Xixia Yingxue (curved) and Hongfeng (erect) of two Paeonia lactiflora Pall. varieties were studied [30]. It was found that there were differences in the degree of secondary cell wall thickening and the number of thickened secondary cell wall layers and 10 key genes related to lignin synthesis were screened by transcriptome sequencing. However, the underlying molecular mechanism of stem variations between Ilex verticillata 'Oosterwijk' and 'Red sprite' is still unclear. In order to explore the mechanism of plant height and stem growth difference between the two varieties, we carried out phenotypic analysis, transcriptome sequencing and endogenous IAA content determination of 'Oosterwijk' and 'Red sprite' in three growth stages. We found that the phenotypes of the two cultivars were quite different. Through comparative transcriptome analysis combined with the results of endogenous IAA content determination, the difference in auxin content and the expression changes of related genes may be the key factors leading to the difference in stem growth between the two cultivars. This can provide a theoretical basis for plant type regulation and the variety improvement of whorled holly and a foundation for the next step of cultivation of whorled holly.

Morphological and Histological Analysis of the Stem in 'Oosterwijk' and 'Red Sprite'
To compare the differences in stem growth between 'Oosterwijk' and 'Red sprite', we analyzed the stem phenotypes in the different growth phases (every 30 days as a harvesting time point) from 29 May 2021 to 31 July 2021 ( Figure 1A). The results showed that there are significant differences in the way of branching and the length of stem between them. Compared with 0 days, the stem diameter of new branches of 'Oosterwijk' increased by 27.5% and 39.7%, respectively, and that of 'Red sprite' increased by 30.3% and 42.7% in 30 and 60 days ( Figure 1B), respectively. In addition, there was no significant difference in stem diameter between the two varieties in the early stage, but in the third period, the stem diameter of 'Oosterwijk' was significantly larger than that of the 'Red sprite'. We also discovered that the length of the new branch of 'Oosterwijk' was longer than that of 'Red sprite' (about 9 cm), but there was no significant difference in the lengths of the new shoots of the two varieties in 0 d vs. 30 d vs. 60 d ( Figure 1C). From subsequent observations, both 'Oosterwijk' and 'Red sprite' appeared to have new shoots in autumn (Supplementary Figure S1). The autumn new shoot of 'Oosterwijk' was greener, thinner and longer than 'Red sprite'.
To further explore the structural differences between the new stems of the two varieties, we observed the cross-section of the new stems through paraffin sections ( Figure 1A). The results showed that compared with the xylem of the 'Red sprite' stem, the xylem of 'Oosterwijk' was thicker and grew faster. At 60 days, the thickness of the xylem of 'Oosterwijk' was 1743 µm, while the thickness of the xylem of 'Red sprite' was 1284 µm ( Figure 1D). We also found that the xylem area of the 'Oosterwijk' stem contains more ducts: about 1.0 times that of the 'Red sprite' ( Figure 1E). The xylem vessel area of 'Oosterwijk' was always higher than that of 'Red sprite' (about 500 µm 2 ) and the area of xylem vessel in 'Oosterwijk' increased significantly at 60 days, but there was no significant difference in 'Red sprite' ( Figure 1F). In addition, the longitudinal section of the stem was observed and it was found that the wood layer of 'Red sprite' was more compact and detailed (Supplementary Figure S2).

Analysis of Differentially Expressed Genes
In order to further explore the mechanism of the stem growth difference between the 'Oosterwijk' and 'Red sprite', 18 gene expression profiles were constructed from the two genotypes with three replicates, and with each at three growth stages (0, 30 and 60 DAG). Each cDNA library produced 31.67-47.36 million Valid Reads and the proportion of Valid Reads reached more than 91%, including 115.38 Gb with an average of 6.41 Gb/library and a total of 30267 transcripts were obtained. The clean reads were compared to the genome of Ilex verticillata for comparative analysis and it was found that the comparison rate between the samples and the reference genome reached more than 88%, that the average comparison rate was 92.68% and that most of the comparison rates reached more than 90% (Table S1).
To understand the changes in gene expression during the growth transition between two cultivars, we identified differentially expressed genes (DEGs) between one time point and the previous time point. According to q-adjust < 0.05, the difference multiple of updown is two, namely |log2 fc| ≥ 1, and the up-and down-regulation of differentially expression genes were obtained (Figure 2A and Supplementary Figure S3). In the stem development phase of 'Oosterwijk' and 'Red sprite', down-regulated genes were more than up-regulated genes. In O_P3 vs. O_P1 and R_P3 vs. R_P1, the number of up-regulated and down-regulated genes had the largest difference, with 1231 and 1108 up-regulated genes, respectively. At the same time, there were 4024 down-regulated genes in 'Red sprite', but only 2838 down-regulated genes in 'Oosterwijk'. However, the number of up-regulated DEGs was more than the number of down-regulated DEGs in O_P1 vs. R_P1, O_P2 vs. R_P2 and O_P3 vs. R_P3. In O_P1 vs. R_P1, 2321 genes were up-regulated while the number of upregulated genes decreased to 1835 in O_P3 vs. R_P3. This appears to be due to fewer downregulated genes at 60 days in 'Red sprite'. The results of the Venn Diagrams illustrated the overlaps and differences between the different comparison groups ( Figure 3B). In O_P1 vs. R_P1, O_P2 vs. R_P2 and O_P3 vs. R_P3, there were 738 common upregulated DEGs and 884, 598 and 677 up-regulated specific DEGs. In addition, there were 423 common down-regulated DEGs and 672, 583, 248 specific down-regulated DEGs ( Figure 2B).

Analysis of Differentially Expressed Genes
In order to further explore the mechanism of the stem growth difference between the 'Oosterwijk' and 'Red sprite', 18 gene expression profiles were constructed from the two genotypes with three replicates, and with each at three growth stages (0, 30 and 60 DAG). Each cDNA library produced 31.67-47.36 million Valid Reads and the proportion of Valid Reads reached more than 91%, including 115.38 Gb with an average of 6.41 Gb/library and a total of 30267 transcripts were obtained. The clean reads were compared to the genome of Ilex verticillata for comparative analysis and it was found that the comparison rate between the samples and the reference genome reached more than 88%, that the average comparison rate was 92.68% and that most of the comparison rates reached more than 90% (Table S1).
To understand the changes in gene expression during the growth transition between two cultivars, we identified differentially expressed genes (DEGs) between one time point and the previous time point. According to q-adjust < 0.05, the difference multiple of updown is two, namely |log2 fc| ≥ 1, and the up-and down-regulation of differentially expression genes were obtained ( Figure 2A and Supplementary Figure S3). In the stem development phase of 'Oosterwijk' and 'Red sprite', down-regulated genes were more than up-regulated genes. In O_P3 vs. O_P1 and R_P3 vs. R_P1, the number of up-regulated and down-regulated genes had the largest difference, with 1231 and 1108 up-regulated genes, respectively. At the same time, there were 4024 down-regulated genes in 'Red sprite', but only 2838 down-regulated genes in 'Oosterwijk'. However, the number of up-regulated DEGs was more than the number of down-regulated DEGs in O_P1 vs. R_P1, O_P2 vs. R_P2 and O_P3 vs. R_P3. In O_P1 vs. R_P1, 2321 genes were up-regulated while the number of upregulated genes decreased to 1835 in O_P3 vs. R_P3. This appears to be due to fewer downregulated genes at 60 days in 'Red sprite'. The results of the Venn Diagrams illustrated the overlaps and differences between the different comparison groups ( Figure  3B). In O_P1 vs. R_P1, O_P2 vs. R_P2 and O_P3 vs. R_P3, there were 738 common upregulated DEGs and 884, 598 and 677 up-regulated specific DEGs. In addition, there were 423 common down-regulated DEGs and 672, 583, 248 specific down-regulated DEGs ( Figure  2B).

GO Classification and KEGG Analysis of Differentially Expressed Genes
To explore the major functional categories and enrichment pathways of DEGs, the GO classification and KEGG analysis of DEGs was performed. The DEGs were annotated in the database. Based on GO analysis, we found that DEGs were annotated to the biological (physiological) process (GO:0008150), regulation of transcription, DNA-templated  Figures S4-S6). There were 280 genes related to auxin or tryptophan among these annotations under BP, CC and MF.

Co-Expression Network Analysis of Genes in the Module
To understand the association between stem difference and phenylpropanoid, phenylalanine more comprehensively and auxin gene regulation, a weighted gene co-expression network analysis (WGCNA) was performed. A total of 24,196 genes were grouped into 25 modules ( Figure 3A). The turquoise module contained the most genes (5959 genes) and the grey module contained the least genes (24 genes) (Table S2). In addition, to find the key modules related to the stem differences of 'Oosterwijk' and 'Red sprite', the correlation between modules and xylem development was analyzed. There were significant differences in xylem development in the third stage, so we chose the modules with the strongest correlation at 60 d. The result indicated that two modules (the. black module and the tan module) were the critical modules and will be identified for further analysis ( Figure 3B).
In the black module, based on GO analysis the majority of DEGs were annotated to protein phosphorylation, defense response and regulation of transcription, DNA-templated under biological process (BP), nucleus, plasma membrane and cytoplasm under cellular component (CC) and protein binding and ATP binding under molecular function (MF) ( Figure 3C). In addition, DEGs were mainly enriched in phenylalanine, tyrosine and tryptophan biosynthesis ( Figure 3E). In the tan module, based on GO analysis, the majority of DEGs were annotated in response to heat, in response to high light intensity and in the regulation of transcription. Furthermore, they were DNA-templated under biological process (BP); the nucleus, cytoplasm and chloroplast were under cellular component (CC) and protein binding and DNA binding were under molecular function (MF) ( Figure 3D). Moreover, DEGs were mainly enriched in protein processing in the endoplasmic reticulum, phenylpropanoid biosynthesis and phenylalanine, tyrosine and tryptophan biosynthesis ( Figure 3F). Therefore, the results of WGCNA suggested that phenylpropanoid biosynthesis and phenylalanine, tyrosine and tryptophan biosynthesis might be the important reason for stem growth differences.

Construction of the Correlation Network of DEGs in the Module
In order to understand the potential regulation mechanism of hub DEGs in two modules, the gene co-expression regulatory networks were constructed according to the data of DEGs from four pathways of phenylpropanoid biosynthesis, phenylalanine, tyrosine and tryptophan biosynthesis, phenylalanine metabolism and tryptophan metabolism. These two networks comprised 21 DEGs and intricate relationships existed among these genes, as shown in Figure 4A,B.
1 in 'Oosterwijk' compared with 'Red sprite' except for gene evm.TU.Chr19.361 in the O_P3 vs. R_P3 group. Therefore, these genes may be the potential core genes for stem growth.
(evm.TU.Chr3.802), sec-independent protein translocase protein like (evm.TU.Chr10.1166), transcription factor bHLH94 (evm.TU.Chr5.1284) and transcription factor bHLH96 (evm.TU.Scaffold62.20, evm.TU.Scaffold51.28) were central to the network ( Figure 4B). There were both positive and negative regulatory relations between them. For example, there was a positive regulatory relationship between evm.TU.Chr4.463 and evm.TU.Scaffold62.20 but there was a negative relationship between evm.TU.Chr4.463 and evm.TU.Scaffold51. 28. In addition, the expression levels of these genes were all lower in 'Oosterwijk' compared with 'Red sprite' except for gene evm.TU.Chr19.361 in the O_P3 vs. R_P3 group. Therefore, these genes may be the potential core genes for stem growth. Accordingly, four hub-screened genes in the black module and seven hub-screened genes in the tan module were selected to construct a general network ( Figure 4C). Correspondingly, three genes including evm.TU.Chr19.361 (HCT), evm.TU.Chr2.706 (AST) and evm.TU.Chr5.1284 (bHLH94) were predicted to be the vital genes of stem growth difference. The bHLH94 belongs to the bHLH transcription factor family, which is the second largest family of transcription factors in plants and participates in multiple developmental Accordingly, four hub-screened genes in the black module and seven hub-screened genes in the tan module were selected to construct a general network ( Figure 4C). Correspondingly, three genes including evm.TU.Chr19.361 (HCT), evm.TU.Chr2.706 (AST) and evm.TU.Chr5.1284 (bHLH94) were predicted to be the vital genes of stem growth difference. The bHLH94 belongs to the bHLH transcription factor family, which is the second largest family of transcription factors in plants and participates in multiple developmental processes of plants. Moreover, HCT is a key enzyme involved in the lignin biosynthesis pathway which affects the biosynthesis of G/S monolignol. Studies have shown that AST is involved in the biosynthesis of phenylalanine, tyrosine and tryptophan and that tryptophan is an important precursor of auxin biosynthesis in plants. There was a positive regulatory relationship between HCT and AST but there was a negative regulatory relationship between HCT and bHLH94. In addition, there was a negative regulatory relationship between AST and bHLH94. In the O_P3 vs. R_P3 group, the expression of AST was higher in 'Oosterwijk' than that in 'Red sprite', but the expression of HCT and bHLH94 were lower in 'Oosterwijk' than in 'Red sprite'. Taken together, these results suggest that bHLH transcription factors, lignin synthesis and tryptophan biosynthesis pathways are key factors in stem growth differences between the two cultivars. Among them, AST may be a key regulatory gene which may affect the synthesis of tryptophan, resulting in the change of auxin content. It is then suggested that the difference in stem growth between the two varieties may also be related to the auxin content.

Analysis of DEGs Related to the Auxin Signaling Pathway and Endogenous Auxin Content Determination
To explore the relationship between stem growth and the expression of auxin-related genes, the expression of the genes involved in the auxin signaling pathway was analyzed. The 29 DEGs related to the auxin signaling pathway were screened and selected for further analysis ( Figure 5). Among these DEGs in 3 phases, the expression levels of 14 DEGs were higher in 'Oosterwijk' at each phase than that compared with 'Red sprite' and 9 genes were more highly expressed in 'Red sprite'. These differential genes are mainly concentrated in auxin synthesis, transport and signal transduction pathways. The SAUR genes are the early auxin response gene, which has the fastest and strongest response to auxin and regulates the synthesis and transport of auxin. Our study showed that the expression levels of the three SAUR genes (evm.TU.Chr7.587, evm.TU.Chr7.6 and evm.TU.Chr10.1539) in 'Red sprite' were higher than those in 'Oosterwijk'. However, the expression levels of GH3 genes (evm.TU.Chr4.2515, evm.TU.Chr4.2519 and evm.TU.Chr4.2514), IAA27 (evm.TU.Chr19.709) and AUX1 (evm.TU.Chr7.218) were higher in Oosterwijk than that in 'Red sprite'. These results suggested that the differential expression of DEGs related to the auxin signaling pathway may affect the variation of endogenous auxin content, which in turn leads to the difference in stem growth between the two cultivars. To further understand the relationship between auxin and stem development, the IAA contents of the stem in both cultivars were determined. The endogenous auxin content was always higher in 'Oosterwijk' than that in 'Red sprite' at each phase ( Figure 6). Thus, these results suggested that the auxin content is positively related to the stem growth in the two cultivars.

Analysis of Auxin Content Determination
To further understand the relationship between auxin and stem development, the IAA contents of the stem in both cultivars were determined. The endogenous auxin content was always higher in 'Oosterwijk' than that in 'Red sprite' at each phase ( Figure 6). Thus, these results suggested that the auxin content is positively related to the stem growth in the two cultivars. Based on the determination of IAA content and transcriptome analysis, the nine DEGs including two GH3.6 genes, GH3, ATRP5/AUX1, IAA27, two SAUR36-like genes, GH3.6-like and AIP 10A5-like were screened and summarized in the auxin signaling pathway ( Figure 5). In 'Oosterwijk', GH3, IAA27 and GH3.6-like were all down-regulated at each phase, while two SAUR36-like genes and AIP 10A5-like were all up-regulated. In 'Red sprite', the up-regulation and down-regulation were basically the same as that of 'Oosterwijk'. In addition, the expression levels of the influx vector ATRP5/AUX1 were higher in 'Oosterwijk' than that in 'Red sprite'. Auxin downstream reaction factors of two GH3.6 genes, GH3 and IAA27, were higher in 'Oosterwijk' than in 'Red sprite', while the expression levels of GH3.6-like, two SAUR36-like genes and AIP 10A5-like were lower in 'Oosterwijk' than that in 'Red sprite'. Therefore, these results sug-

Analysis of Auxin Content Determination
To further understand the relationship between auxin and stem development, the IAA contents of the stem in both cultivars were determined. The endogenous auxin content was always higher in 'Oosterwijk' than that in 'Red sprite' at each phase ( Figure 6). Thus, these results suggested that the auxin content is positively related to the stem growth in the two cultivars. Based on the determination of IAA content and transcriptome analysis, the nine DEGs including two GH3.6 genes, GH3, ATRP5/AUX1, IAA27, two SAUR36-like genes, GH3.6-like and AIP 10A5-like were screened and summarized in the auxin signaling pathway ( Figure 5). In 'Oosterwijk', GH3, IAA27 and GH3.6-like were all down-regulated at each phase, while two SAUR36-like genes and AIP 10A5-like were all up-regulated. In 'Red sprite', the up-regulation and down-regulation were basically the same as that of 'Oosterwijk'. In addition, the expression levels of the influx vector ATRP5/AUX1 were higher in 'Oosterwijk' than that in 'Red sprite'. Auxin downstream reaction factors of two GH3.6 genes, GH3 and IAA27, were higher in 'Oosterwijk' than in 'Red sprite', while the expression levels of GH3.6-like, two SAUR36-like genes and AIP 10A5-like were lower in 'Oosterwijk' than that in 'Red sprite'. Therefore, these results suggested that there were positive correlations between five DEGs (two GH3.6 genes, GH3, ATPR5/AUX1 and IAA27) and stem growth, while there were negative correlations between four DEGs (GH3.6-like, two SAUR36-like genes and AIP 10A5-like) and stem growth.

Water Transport Capacity of the 'Oosterwijk' Stem May Be Stronger Than 'Red Sprite'
Stem length is one of the main factors that affect plant appearance. Further study on the difference of stem growth will be helpful to improve and produce more suitable products for market. In our study, these results indicated that the pith, xylem thickness and vessel area of 'Oosterwijk' were larger than those of 'Red sprite'. Therefore, the mechanism of the stem growth difference finally leads to the phenotypic differences in two cultivars.
A study suggested that the tracheid bridge (a possible pathway between vessels and tracheid to transport water) may be the dominant pathway for water transport between vessels in some species [31]; indeed xylem vessels played an important role in transporting water. In Vitis vinifera L., large vessels were more efficient in the volume of transported water compared to narrower ones [32]. In Rhizophora apiculata, increasing the diameter of the vessels can improve water transport efficiency [33]. In Phyllostachys, the average vessel diameter of Phyllostachys bambusoides Sieb, Phyllostachys glauca f, and Phyllostachys glauca McClure showed a proportional relationship to the water transport efficiency [34]. In our study, the vessel area of 'Oosterwijk' was larger than that of 'Red sprite'. This result suggested that the water transport capacity of 'Oosterwijk' may be stronger than that of 'Red sprite'. In addition, a study showed that late-successional traits V.dulitensis had significantly smaller vessel areas and diameters but higher vessel density [35]. These traits were similar to those of 'Red sprite' in this study. Therefore, these results suggested that 'Red sprite' may be the later successional species and adopted a conservative resource acquisition and growth strategy.

AST, HCT and bHLH94 May Be the Key Genes to the Regulation of Stem Growth
In the O_P3 vs. R_P3 group, the DEGs related to phenylpropanoid biosynthesis, phenylalanine metabolism, the phenylalanine, tyrosine and tryptophan biosynthesis pathway and transcription factors were selected from the black and tan module. In the black module and the tan module, 10 DEGs and 11 DEGs were screened, respectively. Finally, we obtained three hub genes, namely AST, HCT and bHLH94, that might be related to the regulation of stem growth.
In two contrasting maize inbreds viz., I110 (susceptible) and I172 (tolerant), I110 had a higher AST expression level than that of I172, which consumed more energy and oxygen [36]. In Ilex verticillata, 'Oosterwijk' had the higher AST expression level compared with 'Red sprite'. Thus, this result also suggested that 'Oosterwijk', with the high AST expression level, may consume more energy and oxygen. In addition, the study in Dianthus chinensis suggested that the high expression level of HCT may cause high lignin content and huge stems [37]. This finding was consistent with the phenotype of 'Oosterwijk' and its high expression level of HCT. This suggested that HCT may be positively related to stem growth. Moreover, the study showed that the bHLH transcription factor played a role in poplar and tomato stem development [38,39]. In Ilex verticillata, the expression level of bHLH94 in 'Oosterwijk' was lower than that in 'Red sprite', which suggested that bHLH94 may negatively regulate the development of the xylem and stem. These results suggested that AST, HCT and bHLH94, as the center of the network, played an important role in the regulation of stem growth and could be used as candidate genes for the further study of Ilex verticillata stem growth regulation.

Auxin Affected the Growth and Development of the Stem
The study showed that the change in IAA content may be related to stem growth and the auxin concentration directly affected the development of the stem in soybeans [40]. Another study showed that endogenous auxin levels of a novel spontaneous cucumber mutant (named si107) with short internodes were decreased [41]. In Ilex verticillata, the IAA content of 'Oosterwijk' was higher than that compared with 'Red sprite' in three phases. The stem length of 'Oosterwijk' were larger than that of 'Red sprite'. Thus, these results suggested that the content of IAA was positively related to the growth and development of the stem.
A previous study showed that the treatment of white clover with exogenous IAA will result in the down-regulation of the auxin-responsive genes GH3.6 and IAA27 [42]. However, Arabidopsis thaliana treated with NAA induced the expression of GH3.6 [43]. In rice, the study showed that GH3.6 was highly correlated with mesocotyl elongation [44]. In our study, all IAA pathway-related genes were screened for deeply digging out the relationship between the content of IAA and gene expression levels. The result showed that the GH3.6 expression level of 'Oosterwijk' was higher than that of 'Red sprite'. In addition, the study showed that the IAA27 transcript in blueberries was highly accumulated in shoots, but the overexpression of VcIAA27 in Arabidopsis will lead to auxin deficiency and form dwarf plants [45]. Additionally, the study showed that IAA27 in Arabidopsis was a negative regulator of root branching [46]. In Ilex verticillata, the IAA27 content level of 'Oosterwijk' was higher than that of 'Red sprite'. Moreover, some studies demonstrated that SAUR36/RAG1 was related to the regulation of hypocotyl elongation in Arabidopsis thaliana [47]. In our study, the SAUR36 expression level of 'Oosterwijk' was lower than that of 'Red sprite'. Furthermore, in Ilex verticillata, the expression level of ATRP5 was higher in 'Oosterwijk' than that in 'Red sprite' but the expression level of AIP 10A5-like was lower in 'Oosterwijk' than that in 'Red sprite'. These results suggested that GH3.6, IAA27 and ATRP5 may be positively related to the growth and development of the stem, while SAUR36 and AIP 10A5-like may be negatively related to the growth and development of the stem.

Plant Growth Conditions
Ilex verticillata 'Oosterwijk' and 'Red sprite' were cultured at the East-Lake campus of Zhejiang Agricultural and Forestry University in early May 2021. The plant tissues of 'Oosterwijk' and 'Red sprite' were collected for experiments every month from 29 May 2021 to 31 July 2021. Stems close to the bottom of new branches were collected. The stems for the section were collected in the FAA fixative solution, aspirated in the laboratory and stored at 4 • C. The stems for transcriptome and IAA content determination were frozen in liquid nitrogen and stored at −80 • C until the extraction of total RNA. The length and diameter of new branches of Ilex verticillata 'Oosterwijk' and 'Red sprite' were measured and counted at least three biological replicates.

Histological Analysis
The FAA-fixed samples were performed with different concentrations of ethanol and xylene. Then, the samples fixed in FAA were dehydrated and embedded in paraffin. Stems were sectioned for histology using a Vibratome Series 1000 (Heath Company) to a thickness of 50 to 100 µm. Next, dyeing with saffron and fast green was carried out. Finally, they were observed and photographed with a positive fluorescence microscope after sealing the slices. The xylem thickness and vessel area were measured using the Image J software every three slices. The measurement data of xylem thickness came from three randomly selected sections and the vessel area was calculated from three randomly selected sections of a certain size (yellow area in this study).

RNA Isolation and Illumina Sequencing
The total RNA was extracted using Trizol reagent (thermofisher, 15596018) following the manufacturer's procedure. The quantity and purity of total RNA were analyzed by Bioanalyzer 2100 and RNA 6000 Nano LabChip Kit (Agilent, Santa Clara, CA, USA, 5067-1511) with RIN number > 7.0. Then, mRNA was purified from total RNA (5 ug) using Dynabeads Oligo(dT) (Thermo Fisher, Waltham, MA, USA) with two rounds of purification. After purification, the mRNA was fragmented into small pieces using divalent cations under an elevated temperature. Then, the average insert size for the final cDNA library was 300 ± 50 bp. Then, the 2 × 150 bp paired-end sequencing (PE150) was performed on an Illumina Novaseq™ 6000 (LC-Bio Technology CO., Ltd., Hangzhou, China) in accordance with the recommended protocol.

De Novo Assembly and Functional
Firstly, cutadapt software and in-house perl scripts were used to remove the reads that contained undetermined, adaptor contamination and low-quality bases. Then, the sequence quality was verified using FastQC (http://www.bioinformatics.babraham.ac. uk/projects/fastqc/, 0.11.9, accessed on 16 February 2022) including the Q20, Q30 and GC-content of the clean data. All downstream analyses were based on high quality clean data. After that, a total of G bp of cleaned, paired-end reads were produced. Then, we aligned reads of all samples to the Ilex verticillata reference genome using the HISAT2 (https://daehwankimlab.github.io/hisat2/, version: hisat2-2.2.1, accessed on 16 February 2022) package which initially removed a portion of the reads based on quality information accompanying each read and then mapped the reads to the reference genome.

Identification of Differentially Expressed Genes
The differential expression genes (DEGs) between the sample categories were identified using the edgeR package (https://bioconductor.org/packages/release/bioc/html/ edgeR.html, accessed on 16 February 2022). Usually, it is evaluated from two aspects of the difference multiple and the significant level and the differentially expressed genes are screened. Here, we use the difference multiple FC ≥ 2 or FC ≤ 0.5 (|log2 fc| ≥ 1) and p-value < 0.05 as the standard; the genes screened out are considered as differentially expressed genes.

Weighted Gene Co-Expression Network Analysis
The co-expression gene network was constructed using the WGCNA package of software in reference to the tutorial on the WGCNA official website. At first, genes that were not expressed or were poorly expressed in all samples were filtered out. Then, it was determined as to whether there were outliers according to the sample clustering tree diagram. The outliers would be filtered out. According to the filtered data for computing, module information was mined. The correlation between modules, the correlation between modules and samples and the correlation between genes and modules were calculated and then visualized. To further study the gene modules associated with stem differences, the correlation coefficients between module eigengenes and the samples were calculated. Modules with the largest correlation coefficients were selected as target gene modules in the O_P3 vs. R_P3 group. Then, GO analysis and KEGG analysis of DEGs were performed using GO databases and KEGG databases.

IAA Content Determination
The quantification of endogenous IAA was performed according to the method reported previously with some modification [51]. The ground powder of fresh plant tissues (about 200 mg) was extracted for overnight treatment in methanol containing [ 2 H 2 ]-IAA as internal standards. The crude extracts were further purified by the Oasis MCX SPE cartridge which was first activated and equilibrated with methanol, water and 0.05% formic acid (FA). After the samples loading into the cartridge, the cartridge was sequentially washed with 5% FA and water. Finally, IAA was eluted with methanol and analyzed on a liquid chromatography-tandem mass spectrometry system comprising an ACQUITY UPLC (Waters) and Qtrap 6500 system (AB SCIEX) equipped with an electrospray ionization (ESI) source. The separation of IAA was achieved on a Waters ACQUITY UPLC BEH C18 column (100 × 2.1 mm i.d., 1.

Conclusions
In this study, stem morphological observations and comparative transcriptome analysis were performed on stem tissues of 'Oosterwijk' and 'Red sprite' at three reproductive stages to explore the key genes responsible for the differences in stem growth between these two varieties. The O_P1 vs. R_P1, O_P2 vs. R_P2 and O_P3 vs. R_P3 groups contain 3908, 3630 and 2703 DEGs, respectively. The enriched metabolic pathways included phenylpropanoid biosynthesis, phenylalanine metabolism, tryptophan metabolism, plant hormone signal pathway and phenylalanine, tyrosine and tryptophan biosynthesis. The results of endogenous auxin determination showed that the auxin content in the stem tissue of 'Oosterwijk' was higher and the DEGs related to auxin were also identified as candidate genes for stem growth difference. These further illustrated that the change of auxin content and the expression difference of auxin-related genes may be the key factors for the difference between the two varieties of stems. This study can serve as a basis for the identification and functional study of auxin and the stem development of Ilex verticillata and help to illustrate the molecular mechanism of Ilex verticillata stem growth difference, which is expected to breed new Ilex verticillata varieties.