Transcriptional alterations during proliferation and lignification in Phyllostachys nigra cells

Highly-lignified culms of bamboo show distinctive anatomical and mechanical properties compared with the culms of other grass species. A cell culture system for Phyllostachys nigra has enabled investigating the alterations in cellular states associated with secondary cell wall formation during its proliferation and lignification in woody bamboos. To reveal transcriptional changes related to lignification in bamboo, we analyzed transcriptome in P. nigra cells treated with the synthetic auxin 2,4-dichlorophenoxyacetic acid (2,4-D) and the synthetic cytokinin benzylaminopurine (BA) by RNA-seq analysis. We found that some genes putatively involved in cell wall biogenesis and cell division were up-regulated in response to the 2,4-D treatment, and the induction of lignification by the BA treatment was correlated with up-regulation of genes involved in the shikimate pathway. We also found that genes encoding MYB transcription factors (TFs) show correlated expression patterns with those encoding cinnamyl alcohol dehydrogenase (CAD), suggesting that MYB TFs presumably regulate secondary cell wall formation in the bamboo cells. These findings suggest that cytokinin signaling may regulate lignification in P. nigra cells through coordinated transcriptional regulation and metabolic alterations. Our results have also produced a useful resource for better understanding of secondary cell wall formation in bamboo plants.

Bamboo is an ecologically and economically important grass species. It belongs to the largest subfamily, the Bambusoideae, in the grass family (Poaceae) 1,2 , which contains more than 1,500 species that are adapted to diverse climates. It has been exploited for a range of uses such as food, medicine, charcoal, and housing materials, especially in Asia 3 . Owing to their wide utility and productivity, bamboo species are increasingly regarded as a valuable resource for use in renewable energy in the development of a low-carbon society 4,5 .
It is well known that bamboo presents unique biological properties in its vegetative growth and sexual reproduction. It has a rhizome system for lateral growth and forms highly lignified woody culms for longitudinal growth without secondary growth, which are its distinguishing characteristics compared with other grass species and tree species. Moreover, bamboo species often have flowering intervals from several to more than a hundred years, which is another characteristic feature of the sexual reproduction of bamboo species. To elucidate gene regulatory networks involved in these biological phenomena observed in bamboo species, several studies have utilized transcriptome analyses, and identified spatiotemporal expressions of genes explored across different tissues and developmental stages [6][7][8][9] , which improved the understanding of the molecular mechanisms underlying the development and growth in bamboo. However, these analyses provided little information at the cellular level, and did not identify the molecular mechanisms of cellular differentiation associated with its highly-lignified culm formation.
Scientific REPORTS | (2018) 8:11347 | DOI: 10.1038/s41598-018-29645-7 Cell culture systems have been established in some model plant species, such as Arabidopsis T87 10 and tobacco BY-2 11 , and exploited to investigate a wide range of aspects of plant cell biology. Recently, Ogita et al. established a novel xylogenic suspension culture approach in the bamboo Phyllostachys nigra (resource number in RIKEN BioResource Center; rpc00047) that enabled investigation of lignification in living bamboo cells 12 . The cultured P. nigra cells showed cell wall thickening and proliferation in response to treatment with the synthetic auxin 2,4-dichlorophenoxyacetic acid (2,4-D), and lignification occurred in response to treatment with the synthetic cytokinin benzylaminopurine (BA). After 3-5 days of induced lignification, the cells showed xylogenic differentiation, the presence of fiber-like elements with cell wall thickening, and tracheary elements with formation of perforations 12 . Elucidation of the global gene expression profiles of the suspension culture cells under lignification conditions should allow identification of the gene groups important to this process and enable the characterization of gene networks involved in lignification.
The highly conserved genic regions among Phyllostachys species suggest that the draft genome sequence of P. edulis (moso bamboo) 13 can provide a reference genome sequence for RNA-seq-based transcriptome analyses to investigate gene expression patterns in related bamboo species whose whole genome sequences have not yet been deciphered 14 . In-depth analysis of the transcriptome dynamics in response to induced lignification in bamboo cells will provide new insights into the molecular basis of cellular differentiation.
In this study, we aimed to reveal the transcriptional regulatory networks underlying the lignification process of bamboo at the cellular level. We used RNA-seq based transcriptome analysis to obtain an overview of the gene expression of cultured P. nigra cells, rpc00047, and sought to identify the key pathways and transcription factors involved in its lignification process.

Results and Discussion
Overview of the transcriptome analysis of P. nigra cells. We sequenced mRNAs from control and treated P. nigra cells, and found that almost all of the filtered reads could be mapped to the P. edulis draft genome. The P. nigra cells were cultured with treatments of either 2,4-D or BA, and sampled at four and seven days after the initiation of the treatments. Although the cross-platform assessments suggested that Illumina and Ion Torrent would present approximately similar results in RNA-seq based transcriptome profiling, each of them could have platform-specific differentially expressed genes 15 . To minimize biases between the platforms, we applied the Illumina and Ion Torrent sequencing platforms for our RNA-seq analysis of P. nigra cells. From the sequenced mRNAs, we obtained 783 million reads amounting to approximately 78 gigabases in the filtered dataset; 93.22% of these sequences mapped to the P. edulis draft genome (Supplementary Table S1). Thus, even though we used the P. edulis draft genome 13 as the reference sequence, we obtained a high rate of successful mapped reads suggesting that the P. edulis draft genome provides a useful reference genome sequence to analyze transcriptomes in bamboo species, probably due to their conserved genic sequences. We identified 25,443 P. nigra genes significantly expressed in the cells (at least one condition with average RPM values of replicated samples ≥1), which are corresponding to the counterparts annotated in the P. edulis draft genome. These results indicate that, in the P. nigra cells, genes corresponding to as much as 80% of the genes annotated in the P. edulis genome are detectable as significantly expressed genes ( Supplementary Fig. 1a). Comparison of datasets from two duplicate samples after seven days BA treatment and sequenced on the Illumina platform gave Pearson's correlation coefficients (PCC) of up to 0.996. Additionally, comparison of datasets from the same sample conditions using the two sequencing platforms gave high PCC values (e.g., 0.930 between control conditions); the slightly lower PCC values across sequence platforms likely reflect differences in the sequencing methodologies ( Supplementary Fig. S1b). To our knowledge, this is the first study of deep transcriptome analyses of P. nigra, and the data from the study serve as a resource of P. nigra transcripts, which offer clues to identifying genes related to cellular differentiation and lignification in bamboo.

Expression of monolignol pathway genes in response to hormonal treatment of P. nigra cells.
An expression analysis of monolignol pathway genes in P. nigra showed expression of genes putatively involved in the lignification process in the cultured cells. In our previous observation, the P. nigra cells treated with auxins such as 2,4-D or picloram showed increased cell division and suppression of lignification, whereas cells treated with BA showed induced lignification 12 . Moreover, the P. nigra cells under the BA treatment presented increased signals of phloroglucinol-HCl, indicating induction of lignification, and found transcriptional changes in some xylogenesis-related genes including PAL, C4H, CCoAOMT, and CCR induced at day 4 of treatment with BA 12 . To reveal the transcriptional differences underlying the cellular responses against these hormonal treatments observed in the P. nigra cells, we assessed the expression patterns of P. nigra genes putatively involved in monolignol biosynthesis. We found that P. nigra genes encoding CCR and C3H were down-regulated in response to 2,4-D treatment, and that some downstream genes in the monolignol biosynthesis pathway, such as CAD, F5H, and COMT, were up-regulated in response to BA treatment (Fig. 1). Specifically, we found that three genes putatively encoding CAD (homologous to PH01000043G2130, PH01000043G2150, and PH01003504G0010 in P. edulis), F5H (homologous to PH01000012G2270 in P. edulis), and COMT (homologous to PH01000383G0390 in P. edulis) showed a clear response to BA treatment, suggesting their coordinated gene expressions associated with cellular lignification in P. nigra cells. We also found that some genes, such as those encoding PAL, C4H, and 4CL, were up-regulated in response to both 2,4-D and BA treatments. These results suggest that the specific up-regulation of genes encoding CAD, F5H, and COMT in response to BA treatment may presumably be molecular differences associated with the differential cellular responses. For some copies in each gene group, our results are consistent with the cellular responses that initiate differentiation and lignification as well as the expression patterns of genes investigated in the previous study of P. nigra cells 12 . We also found some gene copies, even those encoding the same enzyme that showed different patterns of expression and/or a low level of expression in all conditions, suggesting that subfunctionalization and/or nonfunctionalization may have caused diversification of the expression patterns of these putative paralogous genes. Through our transcriptome analysis, we identified genes involved in the monolignol pathway in P. nigra that were expressed consistently with lignification, suggesting these genes and orthologs will be useful expression markers for monitoring the lignification process in bamboo species. P. nigra gene expression in response to 2,4-D and BA treatments. Comparing the gene expression in cells cultured under control and hormone-treated conditions, we identified a number of DEGs in the latter treatments. We first sought to identify the genes whose expression was responsive to 2,4-D and BA treatments by comparing up-and down-regulated genes. Our results showed that BA treatment triggered a change in expression of a larger number of genes than the 2,4-D treatment (Fig. 2a). Using the threshold adjusted p-value < 1e-03, 1,404 genes were found to be specifically up-regulated in samples from the 4-day treatment with BA compared with that in the control, whereas 177 genes were specifically up-regulated in 4-day treatments with 2,4-D; in addition, 748 genes were up-regulated in both treatment groups (Fig. 2a). In view of clustered gene expression patterns among DEGs, we found that more genes respond significantly to BA than to 2,4-D (Fig. 2b). These results indicate that many genes show altered expression patterns in response to the hormone treatments in P. nigra cells, which suggests that a broad range of cellular systems are influenced by the hormone treatments.

Functional classification of the DEGs. Through enrichment analysis of functional classes and pathways
among the DEGs up-regulated in response to the 2,4-D and BA treatments, we assessed the cellular functions that might be associated with the cellular differentiation induced by the hormone treatments in the P. nigra cultured cells. Among the DEGs specifically up-regulated in response to the 2,4-D treatment, we found that the genes encoding galactose transferase (MapMan #10.3.1.1) and cellulose synthase (MapMan #10.2.1) were enriched in both the 4-day and 7-day treatments (Table 1). We also found that genes putatively encoding cellulose synthase A (KEGG: K10999) were enriched with the 7-day treatment. These findings indicate that P. nigra cells treated with 2,4-D induce expression of the genes related to cell wall biogenesis, which is consistent with our previous observation of thickening and proliferation of the cells in response to the 2,4-D treatment. We also found that laccase activity (KEGG: K05909) was an enriched function among DEGs in response to the 2,4-D treatment. Because secondary wall-associated laccases are required for lignification by catalyzing oxidation of phenolic compounds 16,17 , the 2,4-D treatments may partially activate the process for secondary cell wall formation in P. nigra cells. Among the DEGs specifically up-regulated in response to the BA treatment, we found some enriched functions related to genes encoding enzymes involved in amino acid biosynthesis. Specifically, we found an enrichment of genes involved in the shikimate pathway 18 for biosynthesis of aromatic amino acids (MapMan: #13.1.6.5.1, #13.1.6.5.5, , and cinnamyl alcohol dehydrogenase (CAD) were represented in the monolignol pathway. The expression patterns of P. nigra genes corresponding to their homologs in P. edulis were estimated from the RPM values obtained from the cross-species mapping of P. nigra RNA-seq reads to the P. edulis genome. The color gradient represents normalized gene expression based on z-score of the RPM values. and #13.1.6.1.1, KEGG: K01626), suggesting specific activation of the shikimate pathway in response to the BA treatment ( Fig. 3), which can occur prior to monolignol biosynthesis ( Fig. 1) and the subsequent lignification observed in the P. nigra cells. We also found over-representation of the genes related to transporter activities in the DEGs upregulated in response to the BA treatment (MapMan: #34.99, #34.15, and #34.16, KEGG: K03301 and K03549) ( Table 2), suggesting that BA treatment activates genes encoding transporters and subsequently affects cellular logistics in the P. nigra cells. The list of upregulated DEGs classified to the MapMan binode with the prefix #34 (transport) has showcased genes homologous to various types of transporters, including 12 genes homologous to ATP-binding cassette (ABC) transporters (Supplementary Table S2), which may be involved in the transportation of monolignols 19 . Specifically, four genes encoding putative G family ABC transporters (homologous to PH01000231G0750, PH01002712G0070, PH01002800G0200, and PH01003385G0160 in P. edulis) might be involved in transporting monolignols from the cytoplasm to the cell wall for polymerization in the P. nigra cells. In Arabidopsis, a member of G family ABC transporter, AtABCG29, shows p-coumaryl alcohol transporter activity, and is the first monolignol transporter reported 20,21 . More recently, expression analysis of transporter encoding genes during tracheary element differentiation in cultured Arabidopsis cells suggested that four Arabidopsis ABC transporters; AtABCG11, AtABCG22, AtABCG36, and AtABCG29, may also be involved in lignification as candidate monolignol transporters 22 . The P. nigra cell culture system will provide a useful resource to identify ABC transporters that regulate cellular localization of monolignols in bamboo species, which may offer us novel insights into the evolution of the monolignol biosynthetic pathway in higher plants. In the DEGs upregulated in response to the BA treatment, we also found significant enrichment of a number of genes classified into an unknown functional category (MapMan: #35.2) ( Table 2), suggesting that the BA treatment may affect the expression of genes involved in various cellular functions that remain unexplored. On the whole, these results illuminate the transcriptional alterations of P. nigra cells in response to both the 2,4-D and BA treatments, providing a comprehensive list of genes that may be involved in cellular functions related to proliferation and lignification (Supplementary Tables S3 and S4).
Changes in expression of transcription factor genes in response to BA treatment in P. nigra cells.
We identified transcription factors (TFs) possibly involved in cellular lignification by a comparison of the expression patterns of P. nigra genes putatively encoding transcription factors with those of genes encoding enzymes that catalyze downstream processes in the monolignol pathways, such as CAD, F5H, and COM. This comparison yielded 1,663 genes that putatively encode DNA-binding domains (DBDs) in genes from 60 TF families that are annotated in the P. edulis genome. In addition, based on a comparison of P. edulis and Arabidopsis genomes, we identified genes homologous to Arabidopsis TFs for the promoters involved in cellulose, xylan, and lignin biosynthesis during secondary cell wall formation. These Arabidopsis genes were identified by a yeast one-hybrid assay 23 . Based on the co-expression patterns of these TFs with genes encoding enzymes involved in downstream monolignol pathways, we identified 18 genes putatively encoding TFs, including 7 MYB family genes, three ERF/AP2 family genes, two calmodulin-binding transcription activator (CAMTA), two GRAS family genes, and one bHLH family gene, which showed co-expression patterns with genes putatively encoding CAD, F5H or COMT (PCC ≥ 0.8) ( Table 3). We found that the TF genes were homologous to the AtMYB85 and AtMYB20 genes of Arabidopsis. AtMYB85 is a known lignin-specific transcription factor that regulates lignin biosynthesis genes to activate secondary cell wall formation in Arabidopsis 24,25 . AtMYB20 also regulates secondary cell wall biosynthesis and is induced by NAC transcription factors that regulate secondary cell wall biosynthesis such as SND1, NSTs, and VNDs in Arabidopsis 24 . A co-expression network analysis of genes expressed during internode development in rice identified orthologs of MYB85 and MYB20 as important for secondary cell wall development 26 . These findings suggest that the transcriptional regulatory network for secondary cell wall formation in P. nigra cells might include some TFs conserved between dicot and monocot plants. Moreover, genes in the ERF/AP2 family are homologous to RAP2.12 in Arabidopsis, which is known to have a role in ethylene signaling 27,28 , and possibly regulates the final stages of xylogenesis through ethylene signaling 29,30 . The possible activation of genes for xylogenesis after the induction of lignification in bamboo cells suggests that secondary cell wall formation and subsequent xylogenesis might be coordinated through CK/ethylene crosstalk in bamboo cells 31,32 . We also identified one gene encoding bHLH transcription factors that showed correlated expression with genes for monolignol biosynthesis. This gene was homologous to bHLH105, which encodes IAA-LEUCINE RESISTANT3 (ILR3) that has a crucial role in Fe homeostasis through direct interaction with bHLH34 and bHLH104 33 . It has been reported that both bHLH transcription factors participate in an Arabidopsis gene regulatory network for secondary cell wall biosynthesis 23,34 .

Metabolic differences of P. nigra cells treated with 2,4-D and BA.
In the xylogenic suspension culture, P. nigra cells present differential metabolomic properties during proliferation and lignification in response to treatment with 2,4-D and BA. To reveal metabolic changes and explore its relationship with the transcriptional changes occurring during cellular differentiation, we performed a widely targeted metabolome analysis using CE-MS, with samples of P. nigra cells from the 4-day and 7-day treatments with 2,4-D and BA, respectively, and obtained a metabolome profile dataset composed of accumulation patterns of 214 compounds (Supplementary  Table S5). In the widely targeted metabolome dataset, we found that amino acids synthesized through the shikimate pathway, such as phenylalanine (C_0075: Phe in Supplementary Table S6), tryptophan (C_0102: Trp in Supplementary Table S6), and tyrosine (C_0088: Tyr in Supplementary Table S6), are significantly increased their relative size of MS peaks under the BA conditions (p-value of Welch's t-test < 0.001), suggesting that BA activates the shikimate pathway (Fig. 3), and consequently increases phenylalanine and tyrosine for monolignol biosynthesis in the P. nigra cells (Fig. 2). Moreover, we observed clear metabolic differences between the 4 sample conditions (Fig. 4a,b), suggesting significant metabolic alteration during the cellular differentiation process.
Comparing the metabolome profiles of the P. nigra cells treated with 2,4-D and BA, we identified 120 and 131 metabolites that were differentially accumulated in the cells from the 4-day and 7-day treatments (p-value of Welch's t-test <0.05), respectively, and found that metabolites of amino acids, nucleotide, sugars, and lipids were abundantly accumulated in the cells treated with BA (Supplementary Table S6). Based on the findings of our transcriptome analysis as well as widely-targeted metabolome analysis of the P. nigra cells treated with 2,4-D and BA, we obtained a comprehensive view of the transcriptomic and metabolic alterations occurring in response to the hormonal treatments, which induce proliferation and lignification in a bamboo species (Fig. 5). In response to the 2,4-D treatment, the P. nigra cells activate genes related to cell division and cell growth to promote their proliferation. During this process, they also activate genes associated with biosynthesis of cellulose and hemicellulose, which promote cell wall thickening through primary cell wall formation. In contrast, with BA treatment, P. nigra cells activate genes encoding TFs associated with secondary cell wall formation and the shikimate pathway to synthesize aromatic amino acids, followed by monolignol pathway genes to synthesize monolignol precursors.

Conclusions
Our transcriptome analysis of cultured P. nigra cells that had been induced to undergo proliferation and lignification, identified changes to transcriptional regulatory networks and cellular metabolism, which were presumably related. Functional analyses of the genes encoding TFs that might be involved in lignification in P. nigra will undoubtedly identify regulatory factors for lignification in bamboos. Comprehensive investigation of the lignification process using the P. nigra cell culture system in combination with various -omics analyses will provide a valuable framework for accelerating our understanding of the cellular systems regulating lignification in bamboo species.

Materials and Methods
Cell culture. Bamboo (P. nigra) cells were maintained in suspension culture in modified half-strength Murashige and Skoog (MS) liquid medium 35 supplemented with 3 μM 2,4-D, as described previously 36 . Subcultures were established in 100 ml liquid medium in a 300-ml flask and maintained on a rotary shaker (110 rpm) in the dark at 25 °C. To maintain stable morphology and synchronous growth of the cells, the sedimented cell volume was adjusted to 2.5% every two weeks as described previously 37 . To promote lignification in the cells, 2-week-old cell cultures were transferred to half-strength MS medium supplemented with 10 μM benzylaminopurine (BA) and 3% (w/v) sucrose (lignification conditions) and cultured as described above 12 . Sample preparation for metabolome analysis. The P. nigra cells were immediately frozen in liquid nitrogen and stored at −80 °C until metabolite extraction. Cell samples were weighed and homogenized by Shake Master, BMS-M10N21 (BioMedicalScience, Japan) three times at 1,500 rpm for 2 min, after addition of 500 μl of ice-cold methanol containing 50 μM methionine sulfone as an internal standard. The homogenates were mixed with 500 μl of chloroform and 200 μl of ice-cold Milli-Q water. After centrifugation at 2,300 × g for 5 min at 4 °C, the supernatant was centrifugally filtrated with a Millipore Ultrafree-MC PLHCC HMT Centrifugal Filter Device, 5 kDa (Millipore, Billerica, MA, USA). The filtrate was dried and dissolved in 50 μl of Milli-Q water, and analyzed by CE-TOFMS. RNA extraction. Total RNA was extracted from P. nigra cells using NucleoSpin RNA (Macherey-Nagel, USA), and quality was checked using an Agilent 2100 Bioanalyzer (Agilent, USA). Read processing. The reads from the Ion Torrent-based sequencing that passed the quality control process of the Ion Torrent system were processed by cutadapt 38 to remove sequencing adaptors. The raw reads from the Illumina-based sequencing were trimmed and filtered based on quantity using the FASTX-Toolkit (http://hannonlab.cshl.edu/fastx_toolkit/index.html) with parameter settings of -q 30 -p 80 -v -Q 33.
Reference genome data. The sequence dataset of the draft genome, the coding sequence (CDS), and protein sequences of P. edulis (P. heterocycla var. pubescens) (v1.0) were retrieved from the BambooGDB web site (http://www.bamboogdb.org/) 13,39 . To generate a dataset of structural gene annotations for the P. edulis genome, we mapped the CDS dataset to the draft genome using the GMAP program with default parameter settings, and estimated exon-intron coordinates for the P. edulis genome. A GFF file of the exon-intron coordinates was used to count reads mapped to each gene using featureCount.

IDs of moso bamboo homologs
Closest homologs Promoters Gene expression (z-score of RPM)  Table 3. Transcription factors whose gene expression patterns correlated with the genes involved in monolignol biosynthesis in P. nigra. * Interactions of Arabidopsis TFs for promoters involved in cellulose, xylan, and lignin biogenesis summarized in Kumar et al. 34 based on the Y1H data from Taylor-Teeples et al. 23

Functional annotation and enrichment test.
To predict the functions of P. edulis genes, homology searches were performed using BLASTP (−e = 1e-5, −F = F) against entries of a known protein database (NCBI nr, ftp://ftp.ncbi.nih.gov/blast/db), TIGR Rice Genome Annotation Project (http://rice.plantbiology.msu.edu/), and the protein data present in TAIR release 10 (https://www.arabidopsis.org/). The KAAS web server (http:// www.genome.jp/tools/kaas/) 42 was used to map the protein sequences of P. edulis to metabolic pathways in the KEGG database. In the KEGG pathway mapping, BLAST was used as a search program and "hsa, dme, cel, ath, sce, cho, eco, nme, hpy, rpr, bsu, lla, cac, mge, mtu, ctr, bbu, syn, bth, dra, aae, mja, ape, osa, gmx, and vvi" as search organisms with the single directional best hit method. The Mercator pipeline (http://mapman.gabipd.org/ web/guest) 43 in the MapMan web service was used for functional classification of the protein sequences of P. edulis based on MapMan ontology 44 . The genes putatively encoding transcription factors in P. edulis were annotated based on protein-specific DNA binding domains using Hidden Markov models for 60 transcription factors in  plants with an HMMER search [45][46][47][48][49] . Gene set enrichment analysis (GSEA) of DEGs for MapMan ontology and KEGG pathways was performed by Fisher's exact test.
CE-TOFMS analysis and data processing. CE-TOFMS analysis was performed using an Agilent CE system combined with a TOFMS (Agilent Technologies, Palo Alto, CA, USA) at Human Metabolome Technologies Inc. (HMT, Japan). The samples were diluted two and five times for cation and anion analysis, respectively. Cationic metabolites were separated through a fused silica capillary (50 μm internal diameter × 80 cm length) with Cation Buffer Solution, H3301-1001 (HMT, Japan). Samples were injected at a pressure of 50 mbar for 10 s with the voltage for the CE set at 27 kV. Electrospray ionization-mass spectrometry (ESI-MS) was conducted in positive-ion mode with voltage set at 3 kV. Anionic metabolites were measured through the fused silica capillary (50 μm internal diameter × 80 cm length) with Anion Buffer Solution, H3302-1021 (HMT, Japan). Samples were injected at a pressure of 50 mbar for 25 s with the voltage for the CE set at 30 kV. The ESI-MS was conducted in the negative-ion mode with the voltage set at 3.5 kV. Mass data for the cationic and anionic metabolites were acquired in a range of 50-1,000 m/z. The data were preprocessed using MasterHands software (HMT, Japan). Each metabolite was identified based on m/z and migration time of the MS peak through database search against the HMT database, and was quantified based on the peak area. Differentially accumulated metabolites were identified with a threshold of p < 0.05 in Welch's t-test across the sample conditions. Data Availability. RNA-seq dataset: DDBJ Sequence Read Archive accession number DRA006159.