Transcriptomic analysis reveals transcription factors involved in vascular bundle development and tissue maturation in ginger rhizomes (Zingiber officinale Roscoe)

Ginger (Zingiber officinale Roscoe) is an important vegetable with medicinal value. Rhizome development determines ginger yield and quality. However, little information is available about the molecular features underlying rhizome expansion and maturation. In this study, we investigated anatomy characteristics, lignin accumulation and transcriptome profiles during rhizome development. In young rhizomes, the vascular bundle (VB) was generated with only vessels in it, whereas in matured rhizomes, three to five layers of fibre bundle in the xylem were formed, resulting in VB enlargement. It indicates VB development favouring rhizome swelling. With rhizome matured, the lignin content was remarkably elevated, thus facilitating tissue lignification. To explore the regulators for rhizome development, nine libraries including ginger young rhizomes (GYR), growing rhizomes (GGR), and matured rhizomes (GMR) were established for RNA-Seq, a total of 1264 transcription factors (TFs) were identified. Among them, 35, 116, and 14 differentially expressed TFs were obtained between GYR and GGR, GYR and GMR, and GGR and GMR, respectively. These TFs were further divided into three categories. Among them, three ZobHLHs (homologs of Arabidopsis LHW and AtbHLH096) as well as one DIVARICATA homolog in ginger might play crucial roles in controlling VB development. Four ZoWRKYs and two ZoNACs might be potential regulators associated with rhizome maturation. Three ZoAP2/ERFs and one ZoARF might participate in rhizome development via hormone signalling. This result provides a molecular basis for rhizome expansion and maturation in ginger.


Introduction Introduction Introduction Introduction
Ginger (Zingiber officinale Roscoe), belonging to the genus Zingiber of the family Zingiberaceae, is an important monocot crop and cultivated mainly in China, India and tropical areas. Rhizome, the underground stems of the ginger, is commonly consumed as a spice and medicine in worldwide for its rich in 6-gingerol and essential volatile oils with important medicinal effects. Additionally, fresh ginger rhizome is a famous vegetable in Chongqing and Sichuan, China. Recently, since its vital edible and medicinal values, many researches focused on the rhizome growth and development in ginger Xing et al., 2021). Rhizome development is a complex process, accompanied by rhizome enlargement, maturation, and senescence, which involves a series of changes in cell division and expansion, lignification, metabolite accumulation as well as gene regulation mediated by transcription factors. To date, several studies have revealed variations of bioactive compounds during rhizome growth in varied ginger cultivars (Jiang et al., 2018a;Li et al., 2021;Cheng et al., 2021). However, little information was available about the morphological and physiological features during rhizome expansion and maturation in ginger. Liu et al. (2020) characterized the morphological anatomy of rhizome secretory structures at different developmental stages in ginger. In our previous studies, the lignification caused by the cellulose accumulation and S-lignin increases were illustrated in vertically growing rhizomes (the elongating rhizomes) of ginger (Chen et al., 2020b;Zhang et al., 2022). Rhizome development characteristics are varied dramatically due to different ginger cultivars and different rhizome growing patterns. Because rhizome development correlates with ginger yield and quality, it's of great significance to explore the cytological and physiological events underlying rhizome enlargement and maturation in horizontally growing rhizomes (the swelling rhizomes) of ginger.
Rhizomes are means of asexual reproduction, favoring plants to produce numerous clones effectively. Ginger is a representative species developing rhizome. However, the model species, including Arabidopsis, tomato, and rice, do not develop rhizome, little is known about the molecular mechanisms underlying rhizome initiation and development (Guo et al., 2021). Rhizome growth includes many hormones signaling pathways. GA3 can promote rhizome elongation in tall fescue (Ma et al., 2016). In rhubarb, endogenous ethylene is key regulator of rhizome growth and JA promotes rhizome formation, possibly through inducing endogenous ethylene biosynthesis (Rayirath et al., 2011). In Lotus, hormone signal transduction genes such as GAI, PYL, and ARF are differentially expressed between the elongating rhizomes and the swelling stage rhizomes, indicating crucial roles for GA, ABA, and auxin in rhizome swelling (Yang et al., 2015). In our previous study, we observed that GA, cytokinin and ABA were closely related to cellulose production and rhizome growth in ginger (Chen et al., 2020b). Transcription factors (TFs), also known as trans-acting factors, play important roles in plant growth and developmental processes. According to their conserved DNA binding domains, TFs can be divided into various families with different functions, such as MYB, WRKY, NAC, bHLH, AP2/ERF, TCP, NFY, MADS-box, Homeobox, ARF, GRAS, and bZIP, etc. Several TFs have been identified to be involved in rhizome development. Hu et al. (2011) reported that TFs of YABBY, TCP and homeobox families were specifically expressed in rhizome tips of Oryza longistaminata, indicating potential roles in rhizome development. PpHB1, a homolog of Homeodomain-leucine Zipper (HD-ZIP) transcription factor, was implicated to be responsible for procambial development and rhizome bud formation in bamboo (Wang et al. 2010). In Miscanthus lutarioriparius, 19 TFs belonging to NF-YB, MYB and GRAS families were predominantly expressed in rhizome bud compared with the other tissues, suggesting their involvements in rhizome formation and development (Hu et al., 2017). Yang et al. (2022) analyzed three different rhizome parts (rhizome buds, rhizome nodes, and rhizome internodes) of Poa pratensis var. anceps Gaud. cv. 'Qinghai' (QH) and weak-rhizome wild material (SN) using RNA sequencing, and found that MYB, B3, NAC, BBR-BPC, AP2 MIKC_MADS, BSE1, and C2H2 may be key transcription factors regulating rhizome development, which can interact with multiple functional genes related to carbohydrate, secondary metabolism, and signal transduction, thus ensuring the normal development of the rhizomes. In bamboo 3 rhizomes, most PeGATAs were down-regulated during the rapid-growth of bamboo shoots, and overexpressing PeGATA26 can significantly repress the primary root length and plant height of transgenic Arabidopsis plants . Koo et al. (2013) showed that several TFs, including ethylene response factors (ERFs), AUX/IAA proteins, and MADS-box may play important roles in defining rhizome growth and development in ginger. Xing et al. (2021) obtained a total of 163 AP2/ERFs in Z. officinale, and identified potential rhizome growth regulators ZoERF#064 and ZoERF#006, which were orthologs of AtERF11 and PpERF139, responsible for internode elongation and vessel expansion. We also found several MYBs and NACs associated with lignification and the development of the elongating rhizomes in ginger (Zhang et al., 2022). However, there are no available report illustrating the whole view of transcription factors in the swelling rhizomes of ginger. The molecular features underlying rhizome expansion and maturation in horizontally growing rhizomes of ginger still need to be unraveled.
In this study, to clarify the morphological and physiological characteristics of rhizome enlargement and maturation, we performed an analysis of anatomy as well as lignin content determination in rhizomes at different developmental stages. RNA-sequencing (RNA-Seq) was employed to identify the key transcription factors regulating rhizome growth and development. These results provide important genetic resources for regulating rhizome enlargement and maturation and a theoretical basis for developing new breeding strategies to enhance rhizome yield and quality in ginger.

Plant materials
Ginger (Zingiber officinale Roscoe. cv 'Zhugenjiang') used in this experiment was cultivated in ginger germplasm resource of Chongqing University of Arts and Sciences (Chongqing, China). Ginger rhizomes at three different developmental stages were harvested and illustrated in Figure 1A. Stage 1 rhizomes (GYR) are very young rhizomes, which are newly initiated buds in diameter of approximately 10 mm. Stage 2 rhizomes (GGR) are growing rhizomes, which are middle sized rhizomes in diameter of approximately 20-25 mm. Stage 3 rhizomes (GMR) are matured rhizomes, which are the second order of branching rhizomes in diameter of approximately 40-50 mm. The samples were collected from 12 individual plants and rhizomes from 4 plants were pooled and set as one replicate. Three biological replicates were performed. All samples were placed in tubes, frozen in liquid nitrogen and then stored at -80 °C until further use.

Anatomical structure observation and determination of lignin content in ginger rhizome
The rhizomes from GYR, GGR and GMR were cut into thin slices and stained with phloroglucinolhydrochloric acid. These samples were air dried, grinded into powder, and then the total lignin amounts were analyzed using the method according to our previous study (Chen et al., 2020b). To observe the changes of the anatomical structure in these rhizomes, the samples were cut into small slices, fixed in Carnoy's fluid for 48 h at 4 °C, sectioned (6 μm in thickness), stained with phloroglucinol-hydrochloric acid, and finally photographed by a standard bright field light microscope (Eclipse Ci-L, Nikon Instruments Inc., Tokyo, Japan) with a CCD camera.

RNA-Seq analysis
Total RNA from GYR, GGR and GMR was isolated by using QIAGEN RNeasy Plant Mini Kit according to the manufacturer's instructions. One microgram RNA was used for library construction. Sequencing libraries were generated using Illumina TruSeq™ RNA Sample Prep Kit according to our previous method (Chen et al., 2020b). Briefly, Poly(A) mRNA was enriched from 5 μg total RNA using oligo (dT) magnetic beads and fragmented in a thermomixer compact. Then, the first-and second-strand cDNA were 4 synthesized, end-repaired, A-tailed and ligated with Illumina-specific adaptors. Finally, the cDNA library was constructed by PCR amplification of selected DNA fragments, sequenced on the Illumina HiSeq™ 2000 platform and 150 bp paired-end reads were generated.
The clean RNA-seq data from GYR, GGR and GMR were filtered, processed, de novo assembled by using Trinity (version 2.0.6) (http://trinityrnaseq.github.io/) and TGICL software (version r2013-04-11) (https://sourceforge.net/projects/tgicl/). Subsequently, the contigs and unigenes were subjected to protein sequence database including NCBI non-redundant (Nr), Swiss-Prot, and Pfam databases for alignment and functional annotation by using BLASTx (e-value < 0.00001), respectively. The Pearson's Correlation Coefficient r was used as an indicator to assess the correlation between biological replicates of the samples (Schulze et al., 2012). DESeq2 (Love et al., 2014) was employed to identify differentially expressed genes (DEGs) from transcriptome data. The genes with padj < 0.01 and the absolute value of log2 fold change > 1 were regarded as DEGs.

Screening and classification of transcription factors in ginger
The transcriptome data conducted in this study was used to search for all transcription factors (TFs) in ginger. Firstly, according to the unigenes annotation and TFs classification in the Plant Transcription Factor Database (PlantTFDB, http://planttfdb.gao-lab.org/), the putative TFs in different transcription factor families were screened. Then, TBtools was used to extract the sequences of the TFs (Chen et al., 2020a). ORF Finder and ExPASy Translate tool were employed to identify the open reading frames (ORFs) and amino acids of the putative TFs. The query protein sequences were then subjected to NCBI BlastP analysis. Finally, all candidate TF proteins were confirmed by comparison with the conserved domain against NCBI CDD (https://www.ncbi.nlm.nih.gov/cdd/) and Pfam (http://pfam.xfam.org/). The differentially expressed TFs during ginger rhizome development were also screened out. Moreover, to better and more systematically characterize their functions, the TF families with similar functions are divided into one category.

Bioinformatics analysis of bHLH, WRKY, AP2/ERF transcription factor families
In view of the presence of larger differentially expressed members than other TF families, bHLH, WRKY, AP2/ERF transcription factors were subjected to further bioinformatics analysis. The protein sequences of bHLH, WRKY, and AP2/ERF in Arabidopsis thaliana were downloaded from the PlantTFDB. Clustal W in MEGA X software was employed to perform multiple sequence alignment analysis of Arabidopsis and ginger bHLH, WRKY, AP2/ERF protein sequences by using the default parameters. The conserved protein sequences of these TFs were used for the construction of phylogenetic tree via the Neighbor-Joining method based on a boot-strap test replicated 1,000 times, the p-distance method and pairwise deletion. The online software evolview v2 (https://evolgenius.info/evolview-v2/) was employed to beautify the evolutionary tree. MEME suite (http:// meme-suite. org/ tools/ meme) was applied to predict the conserved motifs of bHLH, WRKY, and AP2/ERF proteins respectively. And the parameters were set as follows: maximum number of motifs, 10; minimum width, 6; and maximum width, 50. Then, the conserved motifs were mapped by using TBtools.

Expression profiling of transcription factors during rhizome growth in ginger
According to their functional classification in the previous studies, the differentially expressed transcription factors were roughly divided into three categories related to hormone regulation, rhizome expansion and growth, and rhizome maturation and aging. The expression profiles of the TFs at three different developmental stages in ginger were obtained from RNA-seq data, and the Tbtools was used to draw heatmaps.

Data availability
The raw transcriptome data were deposited in NCBI Sequence Read Archive under BioProject accession number PRJNA807819.

Rhizome growth and maturation process in ginger
In this study, we observed that significant changes occurred in the microstructure and lignin contents during ginger rhizome development. The ginger rhizome consists of epidermis, cortex, and stele. As the rhizome developmental, the vascular bundle was generated and formed a ring in the pericycle. Many vascular bundles are distributed in the cortex and stele, with a higher density in stele than in the cortex. In young rhizomes, only vessels in vascular bundle were clearly visualized, while in matured ginger rhizomes, tangential and radial cell division appeared, and thus three to five layers of fiber bundle in xylem were formed, resulting in vascular bundle enlargement ( Figure 1B). There were no significant changes in the size of the parenchymatous cells during rhizome growth and expansion. However, these parenchymatous cells arranged loosely in matured rhizomes ( Figure 1B). This indicates that rhizome expansion may be attributed to vascular bundles development in ginger. Sections stained with phloroglucinol-HCl (Wiesner reagent) displayed that the matured rhizome tissues were highly lignified ( Figure 1C). In accordance with this, the result also showed that the total lignin content was remarkably elevated in matured rhizomes than that in newly initiated rhizome buds ( Figure 1D). Error bars represent SE of three biological replicates. Duncan's multiple range test was used to analyse the significance, and the different lower-case letters indicate significant (P < 0.05) differences between samples.

Overview of the transcription factors in ginger rhizomes
To explore candidate transcription factors involved in ginger rhizomes growth and maturation, nine libraries including ginger young rhizomes (GYR), growing rhizomes (GGR), and matured rhizomes (GMR) were established for RNA-Seq and the transcriptome data with high quality was obtained. A total of 79.73 Gb clean reads (~7.72Gb for each library) were generated. After being assembled, 340,939 redundant transcripts were obtained. And then we removed redundancy, 89,172 unigenes were obtained with N50 length of 1528 bp, of which 22.31% were longer than 1000bp. Among them, 31,127 unigenes were annotated against the public protein database. A total of 1,264 transcription factors (TFs) were found in the ginger transcriptomic data, including 24 TF families belonging to MYB, WRKY, AP2/ERF, TCP, bHLH, Zinc finger protein, MADS-box, ARF, AUX/IAA, HSF, bZIP, NAC, Homeobox, NFY, GATA, C3H, Bromodomain, Trihelix, SPL, GRAS, KAN, GRF, YABBY, and Whirly. Several large TF families were presented, such as Zinc finger protein (171) Table S1). The R 2 value between three biological replicates was greater than 0.95 ( Figure S1), indicating strong correlations between samples and high reliability of differentially expressed genes. A total of 35 TFs were found to be differentially expressed between GYR and GGR, among them 3 TFs were upregulated significantly during rhizomes growth, while 32 TFs were dramatically down regulated. 116 TFs were found to be differentially expressed between GYR and GMR, of them, 16 TFs were upregulated significantly during rhizomes maturity, while 100 TFs were dramatically downregulated. 14 TFs were found to be differentially expressed between GGR and GMR ( Figure 2B and 2C). Interestingly, the majority of the differentially expressed TFs were downregulated in response to ginger rhizome development. The results imply that these transcription factors may play a negative regulatory role in rhizome expansion and maturation in ginger.

Identification of transcription factors involved in rhizome growth and expansion
Rhizome size is the major trait that determine the ginger yield, therefore, it's of great significance to characterize the molecular mechanism underlying rhizome enlargement. Herein, the transcription factors associated with rhizome growth and expansion were screened, which includes bHLH, MYB, TCP, GATA, Homeobox, MADS-box, NFY, SPL, Trihelix, GRF, YABBY, C3H, and other Zinc finger proteins. Among them, bHLH is a vital TF family regulating plant cell growth and harbors the most differentially expressed members during rhizomes growth in ginger. In view of this, we focused on this TF family. Based on the RNA-seq data, a total of 69 ZobHLH TFs were observed (Supplemental Table S1). A phylogenetic tree of ZobHLHs and AtbHLH proteins was constructed to investigate their evolutionary relationships (Additional File1). ZobHLHs were divided into 20 subfamilies according to the classifications of their homologs in Arabidopsis ( Figure 3A). Subfamily XV possessed the largest number of bHLH members in ginger (17 genes), whereas subfamilies II, IX and VIIIb harbored the fewest members (only one gene each). Several subfamilies, such as VI, VII(a+b), and VIIIa had no paralogs in Arabidopsis. In many cases, there were different numbers of genes in a given subfamily in Arabidopsis and ginger. Motif analysis showed that three bHLH proteins, including ZobHLH41, ZobHLH56, and ZobHLH60 contained the most motifs (9 motifs) among all ZobHLH proteins, while 11 ZobHLH proteins contained only one conserved motif. It is worth noting that motif1 was exhibited in the majority of ZobHLH genes, considered as the highly conserved bHLH domain ( Figure 3B).
According to the transcriptome data, the expression patterns of genes in these thirteen transcription factor families during rhizomes growth were obtained. The results showed that 17 bHLH genes were differentially expressed during rhizome development, suggesting their involvements in rhizome expansion ( Figure 3C, Supplemental Table S2). Notably, all seventeen differentially expressed bHLHs showed higher expression in GYR and a striking decrease in transcript abundance toward rhizome maturity (from GYR to GMR). In addition, several TFs, which include two MYB genes (MYB86-2 and DIVARICATA), one MADSbox (MADS14), two NFY genes (NFYB3 and NFYA9), four homeoboxes (BEL1-1, HOX9, BEL1-2, and HOX32), two C3H genes (C3H12 and C3H2), and two zinc finger protein (LATE and COL10) were dramatically up-regulated during rhizomes enlargement in ginger. We also observed that two MYB genes (MYB15 and MYB38) and GATA20 displayed significant downregulation and upregulation respectively in GGR compared with those in GYR. C3H33 was upregulated during rhizomes maturation (from GGR to GMR). The rest differentially expressed TFs showed downregulation patterns during rhizomes development and maturation in ginger ( Figure 4). Therefore, most transcription factors might negatively regulate ginger rhizome growth and expansion.  Figure 3. Phylogenetic analysis and expression profiles of differentially expressed bHLHs during rhizome development in ginger. (A) Phylogenetic tree represents the relationships among 183 bHLH proteins of ginger and Arabidopsis. The phylogenetic tree was constructed using MEGA-X software by the Neighbor-Joining method, with a bootstrap test replicated 1000 times, the P-distance and pairwise deletion. The bHLH protein sequences were listed in Additional File1. The different coloured arcs indicate different groups of bHLH domains. The red solid circles represent bHLH domains from ginger. bHLH proteins from ginger with the prefix "Zo" indicate "Zingiber officinale". (B) Conserved motifs distribution of the 69 bHLH proteins in ginger. The conserved motifs of the ZobHLH proteins were identified using MEME suite with zoops (zero or one occurrence per sequence) models of minimum and maximum width of 6 and 50 amino acid residues, respectively and 10 maximum number of motifs. The detailed information for each motif was provided in Supplementary Table S3. (C) Expression profiles of 17 differentially expressed bHLH TFs during rhizome development in ginger. Heatmap was performed by using Tbtools. The first two squares represent mRNA levels in GGR and GMR compared to GYR, respectively, and the following square represents mRNA levels in GMR compared to GGR (from left to right). The expression changes were represented by log2 ratio. Red /blue colors correspond to up-/down-regulation of these genes, and log2 ratio ≥ 1 is considered statistically significant.    Expression patterns of TFs, including MYB, TCP, GATA, MADS-box, SPL, trihelix, homeobox, NFY, GRF, YABBY, Zinc finger protein and C3H during rhizome development in ginger determined by RNA-Seq and illustrated by heatmap. Heatmap was performed by using Tbtools. The first two squares represent mRNA levels in GGR and GMR compared to GYR, respectively, and the following square represents mRNA levels in GMR compared to GGR (from left to right). The expression changes were represented by log2 ratio. Red /blue colors correspond to up-/down-regulation of these genes, and log2 ratio ≥ 1 is considered statistically significant 12

Identification of transcription factors associated with rhizome maturation and aging
Numerous studies have proved that several transcription factors are involved in plant senescence, among which the representative TFs are WRKY, NAC and HSF. Herein, WRKY was one of the largest transcription factor families in ginger. A total of 56 WRKY transcription factors were screened from the ginger transcriptome data, named as ZoWRKY1 to ZoWRKY56 (Supplemental Table S1). A phylogenetic tree of WRKY proteins in Arabidopsis and ginger was constructed to characterize their evolutionary relationships ( Figure 5A). According to the classification of Arabidopsis WRKY proteins, ZoWRKYs were divided into three subfamilies, named Group1, Group2, and Group3. Group2 were then classified into Group2a, Group2b, Group2c, Group2d, and Group2e. The number of WRKY members in each subfamily varied greatly, that the Group2b and Group2c subfamily possessed the largest numbers with 11 WRKY proteins, the Group2e subfamily harbored the least numbers with only 3 WRKY proteins. The MEME motif search tool was employed to identify the conserved motifs of 56 WRKYs ( Figure 5B). Motif1 and motif3 were conserved WRKY DNA binding domain, among which, motif1 was found in majority of the ZoWRKYs (49 out of 56). There were up to seven motifs in three members, including ZoWRKY30, ZoWRKY39, and ZoWRKY56, which were on the same branch in the evolutionary tree. It indicates that WRKY genes in the same group had similar motifs.
The expression profiles of the differentially expressed WRKY TFs from GYR to GMR were shown in Figure 5C. We observed that nine WRKY genes, including ZoWRKY5, ZoWRKY6,ZoWRKY26,ZoWRKY30,ZoWRKY37,ZoWRKY38,ZoWRKY40,ZoWRKY41,and ZoWRKY47, the homologs of WRKY75, WRKY57, WRKY39, WRKY6, WRKY53, WRKY28, WRKY44, WRKY6, and WRKY26 in Arabidopsis respectively, were dramatically downregulated during rhizome maturation ( Figure 5C, Supplemental Table S2). Additionally, nine NAC genes were differentially expressed during rhizome maturation ( Figure S2, Supplemental Table S2). Among them, NAC010 and NAC053, displayed remarkably upregulation in response to rhizome maturation and aging. NAC083 underwent significant down-regulation during rhizome development (from GYR to GGR), but no significant changes were observed as the rhizomes matured (from GGR to GMR). Additionally, the other NAC genes were downregulated responding to rhizome maturity. All differentially expressed HSF TFs showed remarkable increases in their transcripts with the maturation of ginger rhizomes (from GYR to GMR) ( Figure S2). These differentially expressed ZoWRKYs, ZoNACs and ZoHSFs might positively or negatively regulate ginger rhizome maturation and senescence. The phylogenetic tree was constructed using MEGA-X software by the Neighbor-Joining method, with a bootstrap test replicated 1000 times, the P-distance and pairwise deletion. The sequences of WRKY proteins were listed in Additional File2. The arcs with different colors represent different groups (or subgroups) of WRKY proteins. The red pentagram represents WRKY proteins from ginger, WRKY proteins from ginger with the prefix "Zo" indicate "Zingiber officinale". (B) Conserved motifs distribution of the WRKY proteins in ginger. The conserved motifs of the ZoWRKY proteins were identified using MEME suite with zoops (zero or one occurrence per sequence) models of minimum and maximum width of 6 and 50 amino acid residues, respectively and 10 maximum number of motifs. The detailed information for each motif is provided in Supplementary Table S3. (C) Expression patterns of 12 ZoWRKY genes at different stages of ginger rhizome were determined by RNA-Seq and illustrated by heatmap. Heatmap was performed using Tbtools. The first two squares represent mRNA levels in GGR and GMR compared to GYR, respectively, and the following square represents mRNA levels in GMR compared to GGR (from left to right). The expression changes were represented by log2 ratio. Red /blue colors correspond to up-/down-regulation of these genes, and log2 ratio ≥ 1 is considered statistically significant.

Plant hormone signaling-related transcription factors
Plant hormones have a regulatory effect on rhizome growth and development (Chen et al., 2020b). Several transcription factors functioned in plant hormone signal transduction pathway. In this paper, TF families associated with hormone regulation were identified, which includes AP2/ERF, bZIP, ARF, AUX/IAA, GRAS, and KAN. Of which, AP2/ERF was the largest TF family. 69 ZoAP2/ERFs were identified from the transcriptome data (Supplemental Table S1). To better understand their functional significance, a phylogenetic three of AP2/ERF superfamily in ginger and Arabidopsis were constructed ( Figure 6A). According to the taxonomy in Arabidopsis, the 69 ZoAP2/ERF genes were divided into three subfamilies, namely AP2, ERF, and RAV. There were 12 members belonging to the AP2 subfamily named ZoAP2-1 to ZoAP2-12, 54 members belonging to the ERF subfamily named ZoERF1 to ZoERF54, and only 3 members belonging to the RAV subfamily named ZoRAV1-ZoRAV3 (Supplemental Table S1). Motif analysis demonstrated that motif1, motif2, motif3, and motif6 were representative DNA-binding domain in AP2 and ERF proteins. ZoAP2/ERF in the same subclade possessed similar motifs. For AP2 subfamily, ZoAP2-13 and ZoAP2-15 had the most motifs (nine motifs), while ZoAP2-6 and ZoAP2-22 contained the least motifs with only one motif ( Figure 6B).
Based on RNA-Seq data, the expression profiles of the differentially expressed TFs during rhizome development including AP2/ERF, bZIP, ARF, AUX/IAA, GRAS and KAN, were illustrated in heatmap ( Figure 6C and Figure 7). The results showed that there are six ZoAP2/ERF genes, including ZoAP2-3, ZoAP2-9, ZoAP2-10, ZoAP2-12, ZoERF17 and ZoERF34, exhibited dramatic downregulation toward rhizome maturity (from GYR to GMR). ZoERF40 and ZoERF51 showed remarkable upregulation in GMR compared to those in GGR ( Figure 6C). According to Figure 7, two bZIP genes, bZIP44 and bZIP45, underwent significant up-regulation during rhizome development (from GYR to GGR), whereas no significant changes were observed as the rhizome matured (from GGR to GMR). SCL9, a GRAS TF, showed a dramatic upregulation during rhizome maturity (from GGR to GMR), but no obvious changes in its mRNA level during rhizome expansion. SCL27 was strongly upregulated responding to rhizomes maturation (from GYR to GMR). The rest differentially expressed genes belonging to the bZIP, ARF, AUX/IAA, GRAS and KAN TF families were significantly downregulated responding to rhizome enlargement and maturation. In a word, these differentially expressed TFs related to hormone signaling might positively or negatively regulate ginger rhizome development. Figure 6. Figure 6. Figure 6. Figure 6. Bioinformatics analysis and transcript levels of differentially expressed AP2/ERF transcription factors during rhizome development. (A) Phylogenetic tree represents the relationships among 176 AP2/ERF protein of ginger and Arabidopsis. The phylogenetic tree was constructed using MEGA-X software by the Neighbor-Joining method, with a bootstrap test replicated 1000 times, the P-distance and pairwise deletion. The AP2/ERF protein sequences were listed in Additional File 3. The different coloured arcs indicate different groups of AP2/ERF domains. The red solid circles represent AP2/ERF domain from ginger. AP2/ERF proteins from ginger with the prefix "Zo" indicate "Zingiber officinale". (B) Conserved motifs distribution of the 69 AP2/ERF proteins in ginger. The conserved motifs of the ZoAP2/ERF proteins were identified using MEME suite with zoops (zero or one occurrence per sequence) models of minimum and maximum width of 6 and 50 amino acid residues, respectively and 10 maximum number of motifs. The detailed information for each motif is provided in Supplementary Table S3. (C). Expression patterns of 8 ZoAP2/ERF genes in ginger rhizome at different developmental stages determined by RNAseq and illustrated by heatmap. Heatmap was performed using Tbtools. The first two squares represent mRNA levels in GGR and GMR compared to GYR, respectively, and the following square represents mRNA levels in GMR compared to GGR (from left to right). The expression changes were represented by log2 ratio. Red /blue colors correspond to up-/down-regulation of these genes, and log2 ratio ≥ 1 is considered statistically significant.  Figure 7. Expression profiles of transcription factors associated with hormone signaling during rhizome development in ginger. Expression patterns of TFs including bZIP, GRAS, AUX/IAA, ARF, and KAN during rhizome development were determined by RNA-Seq and illustrated by heatmap. Heatmap was performed by using Tbtools. The first two squares represent mRNA levels in GGR and GMR compared to GYR, respectively, and the following square represents mRNA levels in GMR compared to GGR (from left to right). The expression changes were represented by log2 ratio. Red /blue colors correspond to up-/down-regulation of these genes, and log2 ratio ≥ 1 is considered statistically significant

Discussion Discussion Discussion
Ginger is one of the important medicinal and edible plant species in the world. In particular, the ginger rhizome is often consumed as a vegetable in Sichuan and Chongqing of China. In this study, we investigated the rhizome development process of Zingiber officinale cv. 'Zhugenjiang', a famous cultivar in Chongqing. The results demonstrated that rhizome development was accompanied by the enlargement and lignification of rhizome tissues (Figure 1). This study indicates that the appearance, differentiation, and enlargement of vascular bundles are responsible for ginger rhizome expansion, instead of the growth and division of rhizome endodermis cells. This is in accordance with the previous work by Liu et al. (2020) that the stem endodermis of Z. officinale has no meristematic activity, and the vascular bundles are generated from the pericycle, enlarge and push toward the center following with the rhizome development. Herein, the lignin contents were gradually increased as the rhizomes maturation and aging in these horizontally growing rhizomes of ginger ( Figure 1C and 1D). However, in our previous study, there was no obvious change in lignin accumulation during tissue growth and development in the vertically growing rhizomes of ginger (Chen et al., 2020b).  observed that postharvest dehydration stress can significantly enhance lignin production of ginger rhizomes. Therefore, we considered that lignin accumulation in rhizome was determined by the development stage, the type of rhizome growth, and postharvest storage environment. Lignification, as indicated by the enhancement of lignin biosynthesis, is a biomarker for tissue maturation and senescence in ginger rhizomes. Therefore, rhizome development determines ginger yield and quality. It is very important to characterize the molecular mechanism underlying rhizome enlargement and maturation.
Recent studies showed that several important transcription factors have been identified as key elements in coordinating the growth and development, stress responses and secondary metabolism in ginger. Koo et al. (2013) reported that several TFs including one ERF, five AUX/IAAs, two MADS-boxes and one MYB, may play important roles in defining rhizome growth and development of white ginger. Xing et al. (2021) identified the potential ZoAP2/ERF genes regulating rhizome growth, flower development, and abiotic stresses responses. In mango ginger, several TFs, which includes WRKY, MYB, leucine zipper protein, zinc finger and GATA, were identified to be candidate genes associated with resistance to bacterial wilt pathogen (Prasath et al., 2014). There were several TF genes, such as ABF4, WRKY18, WRKY40, and WRKY70, were involved in the defense of ginger to Ralstonia solanacearum infection (Jiang et al., 2018b). In our previous studies, several NACs and MYBs were identified to facilitate ginger rhizome lignification Chen et al., 2020b;Zhang et al., 2022). Cheng et al. (2021) demonstrated that in 'Zhangliang', nine TFs, such as HD-ZIP1, ERF074, GATA4, probably play critical roles in regulating the biosynthesis of 6-gingerol, key secondary metabolite in ginger. However, to date, there were no available information focusing on the comprehensive analysis of transcription factors in Zingiber officinale cv 'Zhugenjiang' and their involvements in rhizome expansion and maturation, particularly the formation and enlargement of vascular bundle. In this study, a complete picture of transcription factors in 'Zhugenjiang' was presented.

Transcription factors regulating vascular bundles formation and rhizome enlargement in ginger
According to previous reports in many plant species, many TF families, such as bHLH, MYB, TCP, GATA, MADS-box, SPL, homeobox, zinc finger protein, were assigned to the category which might regulate vascular bundles development and rhizome expansion. Among them, several TFs displayed dramatic upregulation or downregulation in expanded rhizomes, compared to those in newly initiated rhizomes, suggesting that ginger rhizome growth might be positively and negatively regulated by these TFs. In accordance with the results in Arabidopsis and Oryza sativa (Guo et al., 2005;Gao et al., 2006), bHLH family was one of the largest families of transcription factors in Z. officinale cv 'Zhugenjiang', followed by the bZIP proteins ( Figure 2). However, bHLH is one of the lowest abundant classes in white and yellow ginger (Koo et al., 2013). It indicates that the compositions and numbers of TFs were varied in different ginger cultivars.
In Arabidopsis, bHLH proteins are key regulators for establishment and maintenance of vascular tissue (De Rybel et al., 2013). They found that a TMO5/LHW bHLH heterodimer causes excessive vascular cell divisions, which acts immediately after tissue specification during the very first division of vascular cells in the early embryo, and controls both the establishment of a vascular tissue containing enough cell files and the indeterminacy of this cell population in the growing postembryonic tissue. Another bHLH transcription factor, SAC51-LIKE (SACL), can heterodimerize with LHW, competing with TMO5/LHW interactions, and thus suppress vascular cell proliferation (Vera-Sirera et al., 2015). In addition to controlling cell division, LHW/TMO5-Like1 complex regulates xylem differentiation in the RAM via a novel negative regulatory system, consisting of LHW-T5L1, ACL5, SACL3, and LHW-SACL3, contributes to maintain RAM size and proper root growth (Katayama et al., 2015). Godiard et al. (2011) reported that MtbHLH1 gene, a homolog of AtbHLH096, is involved in controlling nodule vasculature patterning and nutrient exchanges between nodules and roots in Medicago truncatula. In this study, three bHLHs, including ZobHLH51 and ZobHLH67, LHW homologs in Arabidopsis, as well as ZobHLH24, homolog of AtbHLH096, exhibited remarkable declines in their mRNA levels in expanded rhizomes compared to those in newly initiated rhizomes (Supplemental Table S2), indicating their involvements in vascular bundles initiation and rhizome expansion.
Besides bHLH, other transcription factors have also shown to be involved in vascular cell expansion and rhizome development. In tomato, ectopically expressing FSM1, a MYB gene, negatively affects cell expansion, particularly of those cells with the highest potential to expand, such as those residing inner to the vascular bundles in fruit pericarp. They also demonstrated that the FSM1/FSB1/MYBI(DIVARICATA) complex controls differential cell growth in the developing tomato fruit (Machemer et al., 2011). Interestingly, we observed one DIVARICATA gene in ginger displayed dramatic upregulation in expanded rhizomes compared to those in newly initiated rhizomes (Supplemental Table S2), indicating its potential role in rhizome enlargement. Zhong and Ye (2004) suggest that the avb1 gain-of-function mutation of the IFL1/REV gene, belonging to class III HD-ZIP, alters the positional information that determines vascular patterning and organ polarity in Arabidopsis. PeGATA genes were found to regulate rhizome development and bamboo shoot growth partially via GA signaling pathway . In Populus, two MADS-box genes, VCM1 and VCM2, were identified to regulate the proliferation activity of the vascular cambium and secondary growth (Zheng et al., 2021). In Arabidopsis, GRF4 was shown to be essential for both leaf cell proliferation and embryonic development of cotyledons and shoot tip meristems (Kim and Lee, 2006). VvYABBY4 conferred reduced plant stature, dark green leaves, elongated pistil, and reduced size of fruit and seeds (Zhang et al., 2019). Several studies have revealed that TCP genes regulate cell differentiation and cell proliferation. AtTCP15 is primarily expressed in proliferating cells and vascular tissues, and AtTCP15SRDX expression resulted in enlarged trichomes and pavement cells (Li et al., 2012). In short, we considered that the differentially expressed TFs during ginger rhizome growth, including Homeobox-leucine zipper protein, GATA, MADS-box, GRF, YABBY, and TCP, might play vital roles in regulating vascular bundles development. The exploration of the functional significance of these candidate TFs related to rhizome expansion development need to be further analyzed.

Uncover the transcription factors involved in rhizome maturation and senescence in ginger
A large amount of evidence proves that WRKY, NAC, and HSF are involved in the regulation of tissue maturation and senescence. In ginger, tissue lignification and tissue maturation are tightly coupled processes during rhizome development. Teng et al. (2021) demonstrated that CsWRKY13 functioned as a negative regulator for lignin production, as indicated by the reduced lignin content and the transcripts of lignin biosynthesis genes in overexpression Arabidopsis lines. In alfalfa, MsWRKY11 was proved to be a positive regulator for lignin accumulation. Lignin contents were increased in the stem of OE plants, but decreased in SRDX lines (Wen et al., 2021). Hu et al. (2021) revealed that GhWRKY1 was a positive regulator in lignin biosynthesis, which can bind the promoters of lignin biosynthesis related genes GhPAL6 and GhCOMT1, and activate their expressions, thus enhancing total lignin especially S monomers biosynthesis. WRKY transcription factors also play vital roles in fruit ripening. In kiwifruit, AcWRKY40 was responded to ethylene treatment during kiwifruit postharvest ripening. It can enhance the expression of key ripening genes, including AcSAM2, AcACS1, and AcACS2, and bind to their promoters to activate them (Gan et al., 2021). Certainly, numerous reports demonstrated the regulation of WRKY TFs in plant senescence, for example, AtWRKY6, AtWRKY53 and AtWRKY71 act as positive regulators, while AtWRKY70 functions as a negative regulator (Yu et al., 2021). It is worth noting that ZoWRKY37 and AtWRKY70 are belonging to the same branch ( Figure  5A). The expression of this gene was downregulated in matured rhizomes, compared to those in young rhizomes ( Figure 5C), suggesting that ZoWRKY37 might act as negative regulators in rhizomes maturing. Additionally, several WRKYs, such as ZoWRKY6, ZoWRKY26, ZoWRKY30, the homologs of AtWRKY71, AtWRKY11 and AtWRKY6 respectively, also showed differential expressions between GYR and GMR ( Figure  5C), indicating their potential roles in tissue lignification and maturation in ginger rhizomes.
In this study, several NAC and HSF transcription factors were differentially expressed during rhizome development ( Figure S2), suggesting their involvements in response to ginger rhizome maturation. In other species, the roles of these TFs in lignin biosynthesis (Ohtani and Demura, 2019), tissue maturation and senescence have also been extensively presented (Liu et al., 2022). Xu et al. (2015) showed that EjNAC1 was associated with fruit lignification by activating lignin biosynthesis related genes in loquat. Sun et al. (2021) revealed EgNAC141 is a positive regulator of lignin biosynthesis, resulted in stronger lignification, larger xylem, and higher lignin content were observed in overexpression lines in Arabidopsis. In tomato, several NACs, such as Nor-like and SlNAC4, were proved as positive regulators for fruit ripening (Gao et al., 2018). MaNAC42, a transcriptional activator, was involved in the regulation of fruit ripening in banana under oxidative stress. Ectopic overexpression of MaNAC42 in Arabidopsis delays dark-induced senescence in leaves, indicating that MaNAC42 plays a negative role in senescence (Yan et al., 2021). In rice, ONAC022 functions as a stressresponsive NAC with transcriptional activator activity and plays a positive role in drought and salt stress tolerance through modulating an ABA-mediated pathway (Hong et al., 2016). In switchgrass, overexpression of AtLOV1 reduced the leaf angle of the plant by altering the morphology and organization of epidermal cells in the leaf collar region. At the same time, it altered the lignin content and singlet alcohol composition of the cell wall and resulted in delayed flowering time (Xu et al., 2012). In this study, there were two NACs in ginger, NAC022 and NAC035, homologs of OsNAC022 and AtNAC035 respectively, exhibited notable downregulation in matured rhizomes compared to those in young rhizomes (Supplementary Table S2), suggesting that these NACs might function as negative regulators during rhizomes maturing.
Transcription factors associated with hormone signaling are involved in ginger rhizome development Several transcription factors, including AP2/ERF, ARF and AUX/IAA, bZIP, GRAS, and KAN, are involved in the signal transduction pathway of ethylene, auxin, ABA, GA and cytokinin, respectively. Since hormone plays a crucial role in ginger rhizome development (Chen et al., 2020b), it's of great significance to characterize the relationships between these TFs and rhizome enlargement and maturation. Among them, AP2/ERF TFs harbored the most members, which were considered as key regulators for vascular bundle development and tissue maturation. Etchells et al. (2012) demonstrated that plant vascular cell division was maintained by an interaction between PXY and ethylene signaling. Loss of function analysis of ERF109, ERF018 and AtERF1 resulted in plants with inflorescence stems that were characterized by reduced vascular cell numbers, suggesting that these genes promote cell division in vascular meristems. Overexpression of HARDY, an AP2/ERF gene from Arabidopsis, resulted in thicker stems and more xylem rows per vascular bundle in Trifolium alexandrinum (Abogadallah et al., 2011). In tomato, SlAP2a is proved to be a negative regulator of fruit ripening (Chung et al., 2010). Gao et al. (2020) reviewed the regulation of ERFs in fruit ripening, including SlERF6 and MdERF1B. AIL5, an AP2 protein regulated germination and seedling growth in Arabidopsis (Yamagishi et al., 2009). AtERF12 involved in the regulation of seed dormancy downstream of ETR1 and promoting floral meristem identity (Li et al., 2019;Chandler and Werr, 2020). Overexpression of RAP2.6 conferred hypersensitivity to exogenous ABA and abiotic stresses during seed germination and early seedling growth in Arabidopsis (Zhu et al., 2010). In this study, we observed that several AP2/ERFs, including ZoAP2-3 and ZoERF17, the homologs of AIL5 and AtERF12 respectively, displayed significantly differential expression between young and matured rhizomes. RAP2-6, a homolog of RAP2.6, displayed significantly differential expression between growing and matured rhizomes ( Figure 6). We hypothesized that these TFs might control rhizome expansion and maturation. Berleth et al. (2000) pointed out the crucial role of auxin in continuous vascular differentiation and development. Auxin adjusts vascular cell divisions through coordinating activation of SACL and TMO5 (Vera-Sirera et al., 2015). Xu et al. (2019) identified an IAA9-ARF5 module in Populus tomentosa as a key mediator of auxin signaling to control secondary vascular development and xylem formation. In tomato, ARF2 intersects hormonal signals (ethylene, ABA, cytokinin and salicylic acid) to regulate tomato fruit ripening. Silencing of ARF2A resulted in retarded fruit ripening (Breitel et al., 2016). ARF6 and ARF8 in Arabidopsis are required to promote inflorescence stem elongation and late stages of petal, stamen, and gynoecium development (Liu et al., 20 2014). Herein, ARF3 (a homolog of ARF2-like in Musa acuminata) and ARF6 showed extremely higher expression in young rhizomes than those in matured rhizomes (Supplementary Tale S2). It indicates that these two ARFs might play a vital role in regulating rhizome growth and development in ginger. GRAS transcription factors are widely involved in regulating plant development. SlGRAS4 accelerates fruit ripening by regulating ethylene biosynthesis genes and SlMADS1 in tomato . In maize, ZmSCR1 and ZmSCR1h redundantly pattern distinct cell layers in both the leaves and roots (Hughes et al., 2019). In accordance with this, a GRAS gene in ginger, SCR, was highly expressed in newly initiated rhizomes (Supplementary Tale S2), indicating its potential function in rhizome enlargement. Previous studies also showed that bZIPs strongly affect vascular development (Silveira et al., 2007) and tissue maturation (Liang et al., 2020). This study also identified several bZIPs probably correlated with rhizome growth, including three bZIP59 homologs bZIP59-1, bZIP59-2 and bZIP59-3 as well as TGA10.

Conclusions Conclusions Conclusions Conclusions
Herein, we found that in newly initiated rhizome buds, the vascular bundle (VB) appeared with only vessels in it, while three to five layers of fiber bundle in xylem were formed, resulting in VB enlargement with rhizomes maturation. However, no significant changes in parenchymatous cell division and expansion were observed, indicating that VB development is responsible for rhizome expansion. With rhizome matured, the lignin was strongly accumulated, favoring tissue lignification. Based on RNA-Seq data, a total of 35 differentially expressed TFs were obtained between the GYR and GGR groups (3 upregulated and 32 downregulated), 116 differentially expressed TFs were obtained between the GYR and GMR groups (16 upregulated and 100 downregulated), and 14 differentially expressed TFs were obtained between the GGR and GMR groups (6 upregulated and 8 downregulated). These TFs were further divided into three categories. Among them, three ZobHLHs (homologs of Arabidopsis LHW and AtbHLH096) as well as one DIVARICATA homolog in ginger might play crucial roles in controlling VB development. Four ZoWRKYs and two ZoNACs might be potential regulators associated with rhizome maturation. Two ZoERFs and two ZoARFs might participate in rhizome development via hormone signaling. Their roles would be characterized via genetic and biochemical approaches in vitro and in vivo in our future work.

Authors' Contributions Authors' Contributions Authors' Contributions Authors' Contributions
Conceptualization, Z C, F X and N T; Data curation, Y L and Z L; Formal analysis, Z C and J Z; Funding acquisition, N T; Investigation, Z C and X S; Methodology, Z C and P W; Resources, P W; Software, J Z and L Z; Supervision, J Y and Z C; Visualization, P W and Y H; Writing -original draft, Z C; Writingreview & editing, F X, Z C and N T.
All authors read and approved the final manuscript.
Ethical approval Ethical approval Ethical approval Ethical approval (for researches involving animals or humans) Not applicable.