Determining the different mechanisms used by Pseudomonas strains to cope with minimal inhibitory concentrations of zinc via comparative transcriptomic and metabolic analyses

Background Pseudomonas is one of the most diverse bacterial genera identified in soil, water, plants and clinical samples. Genome sequence analysis has indicated that this genus can be clustered into three lineages and ten groups. Each group can adopt different mechanisms to thrive under zinc-depleted or high-zinc conditions, two environments that are frequently encountered during their environmental propagation. Results The response mechanisms of three prominent Pseudomonas strains ( Pseudomonas aeruginosa PAO1, Pseudomonas putida KT2440, and Pseudomonas fluorescence ATCC13525T) under high-zinc condition were compared using RNA-seq and ultra-performance liquid chromatography–tandem mass spectrometry analysis. Results demonstrated that the three strains shared only minimal similarity at the transcriptional level. Only four genes responsible for zinc efflux were commonly upregulated. P. aeruginosa PAO1 specifically downregulated the operons involved in siderophore synthesis and the genes that encode ribosomal protein, while upregulated the genes associated with antibiotic efflux and cell envelope biosynthesis. The membrane transporters in P. putida KT2440 were globally downregulated, indicating changes in cell permeability. Compared with P. aeruginosa PAO1 and P. putida KT2440, the most remarkable transcriptional variation in P. fluorescence ATCC13525 is the significant downregulation of the type VI secretion system. Metabolite quantitative analysis showed that low concentrations of the metabolites involved in central carbon metabolism and amino acid synthesis were detected in the three strains. Dipeptides containing branched-chain amino acids seems played an important role in alleviating stress caused by zinc ions.

3 three strains apparently used different pathways to reduce zinc entry. In addition, zinc treatment can increase the difficulties of scavenging P. aeruginosa from its colonisation area, and reduce the effectiveness of P. fluorescence as a biocontrol agent.

Background
Anthropogenic industries and domestic activities directly or indirectly release vast amounts of metal ions into the environment. The regional accumulation of these metal ions exerts various effects on ecological systems and human health. Many of these metal ions are toxic metals, such as Al, Bi, Cd, Cr, Pb, Sn, and Ti, which have no known positive biological role [1,2]. Although some transitional metal ions are essential nutrients in all three kingdoms of life, excessive amounts of these ions can also exhibit toxicity [1]. Zn 2+ is the second most abundant essential metal ion after Fe 2+ /3+ , which directly interacts with RNA polymerase, zinc-finger proteins and superoxide dismutase to maintain their functionality [3,4]. DNA-associated zinc potentially influences gene expression by modifying DNA structure and stability [3,5]. Free Zn 2+ in the cytoplasm is normally maintained at the femtomolar to picomolar level [6]. Once this threshold is surpassed, excess Zn 2+ can replace other metal ions from proteins and alter the stability of biomolecules [7,8]. Although zinc is a redox-inactive metal under normal conditions, high levels of intracellular zinc clearly induce the generation of reactive oxygen species (ROS) as evidenced by the increased amount of oxidised protein [9].
Prokaryotes are single cell organisms that have direct contact with their surrounding environment and are separated from the environment only by a double layer cellular membrane. Thus, these organisms are extremely vulnerable to environmental stresses [10]. In response to metal stress, many bacteria have evolved numerous mechanisms, including mobilisation, precipitation, biosorption and biogeochemical cycling, to respond 4 to the adverse effects of metal ions [2]. Pseudomonas is a ubiquitous bacterial genus with high adaptability to various environments [11]. Many Pseudomonas species exhibit high capability to maintain the equilibrium of cytoplasmic metal ions, which directly or indirectly affects metal availability and distribution [12,13]. In our previous study, the transcriptional responses of P. putida KT2440 to low, intermediate, and high levels of extracellular zinc stresses were thoroughly analysed [7]. The results showed that genes associated with metal transport and membrane homeostasis were significantly induced at the lowest zinc level. As zinc concentrations increased, the respiration chains, amino acids and carbon source metabolisms, as well as the biogenesis of Fe-S clusters were adjusted.
However, Pseudomonas comprises approximately 200 species [14], which exhibit considerable diversity in their ecological habitat and physiology [2,15]. Recent phylogenetic analysis has clustered Pseudomonas into three lineages and ten groups [15]. Each Pseudomonas group probably uses different way to alleviate metal-induced intracellular disorder. A comparision of our transcriptome data with the results obtained by Alhasawi et al. [9] showed partial differences between the response of P. putida and P. fluorescence to zinc. For example, the expression of NADPH-generating enzymes, including isocitrate dehydrogenase-NADP + , malic enzyme (ME), and glucose-6-phosphate dehydrogenase (G6PDH), was markedly increased in P. fluorescens. Meanwhile, the most upregulated NADPH-generating enzyme in P. putida KT2440 was ferredoxin-NADP + reductase and the transcription of ME and G6PDH remained unchanged. In addition to P. fluorescence and P. putida, P. aeruginosa represents another important phylogenetic group in the Pseudomonas genus [15,16]. As a widespread environmental bacterium and an opportunistic animal pathogen, P. aeruginosa also encounters changing levels of zinc concentrations [17]. To date, studies on zinc homeostasis in P. aeruginosa mainly concentrated on zinc efflux and zinc associated virulence or antibiotic resistance [3,17,18], other cellular responses remain poorly understood.
In the current study, three prominent Pseudomonas strains (P. putida KT2440, P. aeruginosa PAO1, and P. fluorescence ATCC13525 T ) were challenged with external zinc at minimal inhibitory concentrations (MICs). Their responses were compared at the transcriptomic and metabolic levels based on RNA-seq and ultra-performance liquid chromatography-tandem mass spectrometry (UPLC-MS/MS) to provide a comprehensive understanding of the genes, the cellular processes and the cellular molecules that play crucial roles in zinc adaptation.

MICs of zinc sulphate
Given that metal ions exert does-dependent effects on bacterial cells, comparing the cellular responses of different strains at the same zinc stress level is important. In this study, the MIC of ZnSO 4 was used in the treatment. P. putida KT2440 was highly tolerant to zinc and exhibited a zinc MIC value of up to 1.1 mmol L − 1 . P. fluorescence was the least zinc-tolerant strain with an MIC value of approximately 0.4 mmol L − 1 . P. aeruginosa PAO1 had intermediate levels of zinc tolerance. A cation-defined medium (CDM) plate that contained more than 0.7 mmol L − 1 zinc completely inhibited its colony formation. Therefore, in the following RNA-seq and UPLC-MS/MS analysis, 0.7, 1.1, and 0.4 mmol L − 1 zinc sulphate were used to compare the zinc-induced cellular responses in P. aeruginosa PAO1, P. putida KT2440, and P. fluorescence ATCC13525.

Global transcription features of the three Pseudomonas strains
The numbers of differentially expressed genes (DEGs) during the zinc treatment of P. aeruginosa PAO1, P. putida KT2440, and P. fluorescence ATCC13525 T are provided in Table 1. The comparison of the zinc-treated P. aeruginosa PAO1 with their control samples identified 529 genes that were differently transcribed. Amongst which, 406 genes, 7.1% of the total coding DNA sequence (CDS) in the genome, were upregulated, whereas 123 (2.2%) genes were downregulated. The comparison of the zinc-treated P. putida KT2440 with their control bioreplicates showed that, the expression of 607 genes was significantly altered, with the majority (344 genes, 6.0%) displaying significant induction in transcription compared with those (263 genes, 4.6%) demonstrating decreased transcription. The lowest number of DEGs was observed in P. fluorescence ATCC13525, with only 48 induced genes and 19 repressed genes. In general, KEGG analysis showed that large quantities of DEGs identified in P. aeruginosa PAO1 and P. putida KT2440 can be clustered into the metabolism class ( Fig. 1a), particularly the carbohydrate and energy metabolism subclasses. As an opportunistic animal pathogen, P. aeruginosa PAO1 also regulated many genes belonging to the human disease class. The cellular response of P. fluorescence ATCC13525 differed from those of P. aeruginosa and P. putida, in which approximately 50% of DEGs with clear KEGG classifications were categorised into the environmental information processing group. To better compared the cellular responses occurred in the three Pseudomonas strains, the DNA sequence of the DEGs was extracted from each genome, and then blast against the DEGs identified in the other two genomes with a match cutoff of 70% and an E-value exponent cutoff of 1-e5. As shown in Fig. 1b, only four orthologues were commonly regulated by the three Pseudomonas strains, and 1044 genes were strain-specifically regulated. The transcriptional patterns of P. aeruginosa and P. putida presented high similarities with 58 commonly regulated DEGs. In accordance with the KEGG analysis, P. fluorescence stands out in in all pair-wise blast analysis. Only 19 orthologues were commonly regulated by P. fluorescence and P. aeruginosa, and this value was reduced to 7 in P. fluorescence versus P. putida analysis. Table 1 The quantity of DEGs in P. aeruginosa PAO1, P. putida KT2440, and P. fluorescence ATCC13525 T . Functional analysis of DEGs identified in more than one strain All the four DEGs (cadA, cadR, czcR, and czcS) commonly upregulated by the three strains were involved in metal efflux. cadA has been recognised as an efficient P-type ATPase that moves several divalent heavy metal ions out of the cellular membrane [19], whilst czcRS has been proven as an important two-component system that connects metal efflux and entrance tunnels. In P. aeruginosa, unphosphorylated CzcR decreased the expression of outer membrane porins, OprD. By contrast, the presence of phosphorylated CzcR induced the most effective metal efflux system, czcCBA [ 20]. DEGs mutually regulated by P.
aeruginosa PAO1 and P. putida KT2440 are listed in Additional file 1 (Table S1). As expected, czcCBA operons, which are known to be induced by several divalent metal ions, were found to be highly responsive. The analysis of other DEGs with known functions showed that most of these DEGs can be manually clustered into three categories. The genes involved in membrane structure and channels constituted the largest functional category. In this group, seven genes, including desA, ompB, dgkA, warB, warA, amgR and PA4819 (PP_0034), were strongly upregulated (> 5 fold), whereas oprD and kguT were mutually downregulated. In the other two groups, genes associated with central metabolism and those responsible for protein folding and degradation were generally upregulated. For P. fluorescens, the most prominent response mutually identified in P.
fluorescens and P. aeruginosa was the upregulation of acyl-CoA dehydrogenase (RS07420, RS06160) and protein-tyrosine-phosphatase (RS10450). A ton-dependent receptor (RS23050) and a non-ribosomal peptide synthatase were markedly down-regulated 8 (RS27315) (Additional file 2: Table S2). In addition to the four DEGs commonly regulated by the three strains, P. fluorescens and P. putida shared three more transcripts (Additional file 3: Table S3). Two of them were responsible for active transmembrane transport (RS21740, RS19005). The induction of asparaginase (RS29965) was also identified. The genome locations of these genes that are regulated by more than one strain are spread throughout the genomes, no obvious genome island was observed (Fig. 1c).
DEGs specifically identified in P. aeruginosa PAO1 The DEGs specifically identified in P. aeruginosa PAO1 are shown in Additional file 4 ( Table   S4). The fold-change data indicated that more energy were required for P. aeruginosa to alleviate zinc-induced toxicity. The proton-translocating NADH:ubiquinone oxidoreductase complex (nuoGHIJLMN) and the pyruvate dehydrogenase complex (PA3416-3417), which transform pyruvate into acetyl-CoA were upregulated by approximately 2.3-fold, whilst the synthesis of ribosomal protein was significantly downregulated (rpsB, rpsS, rpsT, rpsR, rpsF, rplS, rplW, rplJ, and rplU). A distinct group of operons involved in lipid A synthesis/modification was highly induced (arnBCADEF), suggesting that the P. aeruginosa cells were in a cell-envelope-stressed state during growth. In support of this, there was also significant upregulation of other genes previously linked to bacterial cell envelop homeostasis, such as glmM (phosphoglucosamine mutase) and pagL (lipid A 3-Odeacylase). The direct raplacement of iron from their binding sites has been recognised as an important mechanism for zinc to exhibit cytotoxicity. Therefore, it is reasonable that the biosynthesis of pyoverdine was downregulated (pvdQAPNOFEDJ) to reduce the unexpected intrusion of zinc ions.
Given that P. aeruginosa is an opportunistic pathogen, monitoring the transcriptional variation of pathogenicity related genes is important. The MIC of zinc strongly enhanced the antibiotic resistance of P. aeruginosa. Two important multidrug efflux systems, mexRAB-oprM [ 21] and mexXY [ 22], were strongly upregulated. Moreover, one operon involved in alginate synthesis (algU-mucABC) was significantly activated. The overproduction of exopolysaccharide alginate caused mucoid conversion in P. aeruginosa [23,24], which increased bacterial metal tolerance via metal chelation [25]. Amongst all the DEGs specifically identified in P. aeruginosa, the most upregulated gene was ptrA. The function of PtrA was studied by two groups, but contrasting conclusions were reported.
One group showed that PtrA directly binds to ExsA, which in turn, suppresses the expression of the type III secretion system (T3SS) [26]. By contrast, the other group demonstrated that PtrA is a periplasmic protein, the expression of which increases the Cu tolerance of P. aeruginosa without affecting basal ExsA [27]. Our data supported that PtrA is not a T3SS repressor because the transcription of the T3SS gene clusters remained unchanged when ptrA was considerably induced by zinc.
DEGs specifically identified in P. putida KT2440 In our previous study, the transcriptional response of P. putida KT2440 to stress-inducing concentrations of zinc was analysed. The results showed that different zinc stress levels strongly influenced (> 4 fold change) the transcription of genes from four functional groups, including metal transporting genes, genes associated with membrane homeostasis, antioxidant-encoding genes and genes involved in basic cellular metabolism [7]. In this study, a twofold change was used as a criterion, and a total of 545 genes were identified as DEGs specifically regulated in P. putida KT2440. After analysing the function of differentially expressed operons, the results further showed that most of the nutrient uptake transporters (about 64% of the total downregulated operons) were downregulated (Additional 5: Table S5), indicating that P. putida KT2440 tended to decrease the permeability of the cell envelope. Only the transporters responsible for methoine, sulphate and sulfonates uptake were upregulated, suggesting that, unlike P. aeruginosa PAO1, P. putida KT2440 did not suffer from a severe iron metabolism perturbation, but the sulphur-containing molecules were considerably disrupted. The synthesis of ribosomal protein was another remarkable difference between P. aeruginosa PAO1 and P. putida KT2440. In P. aeruginosa PAO1, the transcription of ribosomal protein encoding genes was generally downregulated. Meanwhile, their orthologues in P. putida KT2440 were highly stable, and even the transcription of rpmH was increased. In addition, the MIC of zinc significantly increased the transcription of the nickel (nikABCDE) and arsenic (arsBR) efflux systems in P. putida KT2440. Such phenomenon was not observed in the other two strains.
DEGs specifically identified in P. fluorescens ATCC13525 Compared with the transcriptional responses of P. aeruginosa PAO1 and P. putida KT2440, the downregulation of type VI secretion system (T6SS) was the most remarkable feature of P. fluorescens ATCC13525 under zinc stresses (Additional file 6: Table S6). T6SS is widely distributed across diverse bacterial species; around one-third of the sequenced Gramnegative bacteria possess T6SS-associated genes [28]. Bacterial T6SS functions as a contractile nanomachine that delivers effectors upon direct contact with a target cell [29].
Given that these effectors have different functions but frequently disturb the cellular structure, such as the cell wall, nucleic acid, or membrane compartment, T6SS is currently perceived as an antibacterial weapon [28]. In recent study, the T6SS4 in Burkholderia thailandenss was found to secrete a proteinaceous zincophore, which interacts with the outer membrane heme transporter to import zinc under oxidative stress [30]. Given that the transcriptional changes of other zinc importers, such as znu and zip, were not observed in the current study, we postulated that T6SS also plays an important role in zinc transport in P. fluorescens, the downregulation of which reduced the transport of zinc across the outer membrane.
Reverse transcription quantitative real-time polymerase chain reaction (RT-qPCR) validation To confirm the RNA-seq data, six genes were selected from each strain, and their transcription was analyzed via RT-qPCR with three biological replicates. As shown in Fig. 2, the dynamic transcription patterns of all the genes were consistent with the RNAseq data (Fig. 2a), and the observed fold changes for each gene were moderately correlated (R 2 ranged from 80.1-89.6%) (Fig. 2b). Therefore, the RT-qPCR results validated the accuracy of the RNA-seq data. Moreover, previous studies on the transcriptomic profiles of P. putida KT2440 were useful in further confirming the results obtained in this study.

Metabolite changes in the three Pseudomonas strains
Considering that large quantities of DEGs identified in the transcriptome analysis were genes associated with metabolism, we detected the concentrations of metabolites in these bacterial cells. Metabolic profiles of whole-cell extracts were obtained from zinc-treated and control samples exponentially cultured in CDM. Six replicates were performed for each sample. A total of 83 annotated metabolites were significantly changed, and approximately 76% of which were down-regulated. Only two metabolites were commonly regulated by the three Pseudomonas strains (pyroglutamic, 2-phosphoglycerate). The metabolite profiles of P. aeruginosa and P. putida shared more similarity, 18 metabolites were commonly regulated (Fig. 3) putida were only about 70% of that calculated in the cells without zinc treatment, and the transcription of eda (2-keto-3-deoxy-6-phosphogluconate aldolase) were downregulated about 1.7 and 1.5-fold in P. putida and P. fluorescens, respectively. Central carbon metabolism is the main pathway that satisfies the energy requirement of cells. Moreover, these metabolic processes provide important precursors for primary and secondary metabolites. Therefore, the metabolites of many other pathways were also reduced. For example, a wide range of amino acids (ornithine, Lys, Pro, Asp, Trp, and His) or intermediates involved in amino acid biosynthesis were decreased in P. putida. And lower concentrations of metabolites involved in Arg, Glu, Lys, Cys, and Met metabolisms were observed in P. aeruginosa.

Discussion
The genus Pseudomonas is one of the most diverse bacterial genera, comprising approximately 200 species isolated from sources ranging from vegetation to water and contaminated soil to animal clinical samples [2]. With relatively large genomes (typically between 6 Mbp and 7 Mbp), many of these species have evolved delicate gene regulation systems that allow them to thrive under various environmental stresses [16,31]. By combining the transcriptomic and metabolic data obtained in the present study, an overview of the global effects of zinc on three representative Pseudomonas strains was provided.
13 P. aeruginosa PAO1 was originally isolated from a human wound and widely studied as a model opportunistic human pathogen [32]. Soil-originating P. putida KT2440 is the bestcharacterised strain from this species. It is used worldwide as the primary strain for genetic and physiological studies [33]. High genetic diversity was observed among bacteria classified as P. fluorescens, which can be further clustered into three smaller taxonomic subclades [34]. Strains classified into one subclade shared a higher amount of conserved genes than the species-complex as a whole. Comparative genomic analysis showed a high level of genome similarity between P. putida KT2440 and P. aeruginosa PAO1, about 85% of the CDS are shared [16]. Although P. fluorescens ATCC13525 behaved as an out-group amongst the three strains, our bioinformatic analysis indicated that there is still 29% CDS are shared between P. fluorescens ATCC 13525 and P. putida KT2440, and 33% CDS are shared between P. fluorescens ATCC 13525 and P. aeruginosa PAO1. By contrast, the transcriptomic features of the three strains in response to the MICs of zinc exhibited significant divergence. Only 4 orthologues were commonly regulated, and all these orthologues were involved in metal efflux. The exclusion of cytoplasmic ions by membrane efflux pumps has been recognised as an important strategy for prokaryotic cells to withstand metal toxicity [35]. However, this process is not energy-efficient, and damages can happen before these efflux systems are fully expressed. Reducing the entry of excess metal ions can more effectively prevent the generation of metal-induced cytotoxicity. Based on our data, the three strains use different way to preclude the entry of excess zinc ions. P. aeruginosa PAO1 significantly decreased the secretion of siderophore, which was also found to complex with zinc, nickel, cobalt and magnesium as well as iron [36,37]. P. putida universally reduced the expression of membrane transporters. Cell permeability can be significantly decreased. P. fluorescence ATCC13525 has the lowest MIC value of zinc, and a relatively small number of genes were 14 transcriptionally changed. Among these genes, T6SS gene cluster was remarkably downregulated, whose analogues in B. thailandenss play an important role in zinc uptake under oxidative stresses [30]. In addition to the change in zinc transport, it interesting to note the side effect caused by zinc ions. The increased alginate synthesis can facilitate bacterial adherence, biofilm formation, neutralize oxygen radicals, and act as a barrier to phagocytosis [38]. Coupled with the upregulated antibiotic efflux systems, the scavenging of zinc-exposed P. aeruginosa can be more difficult. P. fluorescens has been studied extensively as a plant growth promoter that synthesises toxic metabolites against phytopathogenic microorganisms, and enhancing nutrient availability in soil [39]. The downregulated transcription of T6SS can reduce their competence to occupy phyllosphere and rhizosphere, which in turn, decreases the effectiveness of P. fluorescens to be used as a biocontrol agent.
As transcriptome analysis showed that large quantities of DEGs identified in P. aeruginosa, P. putida, and P. fluorescens can be categorised into the metabolism KEGG group, Hence, different metabolic intermediates were measured in detail using untargeted metabolomics analysis to further explore the mechanisms underlying the effects of zinc exposure on the metabolic pathway. Two classes of metabolites were significantly decreased according to UPLC-MS/MS analysis, including central carbon metabolism and amino acids synthesis, which may be the direct reason why the growth of these bacterial cells stopped. Dipeptide is another group of metabolites that identified at lower concentrations after zinc treatment, especially in P. aeruginosa (Pro-Arg, and Ile-Leu) and P. putida (Pro-Val, Val-Ala, and Ile-Ala). High percentages of branched amino acids were observed, indicating that these dipeptides were not generated by protein degradation. The changes of Ile/Leu/Val containing short peptides were also observed by Jousse et al. [40]. Their metabolomics data revealed that the concentrations of Ala-Leu, Leu-Leu, Ala-Arg, Pro-Phe, Pro-Ile, Leu-15 Phe, and Val-Val were significantly decreased in cold-shocked P. syringae. Dipeptides containing a branched-chain amino acid have been proven to exhibit higher antioxidant activities than dipeptides constituted by other amino acids [41]. Given that high concentrations of zinc ions or low temperature clearly increased the ROS levels in prokaryote cells [42], we inferred that a large amount of branched-chain amino acid containing dipeptides were consumed by ROS quenching. Therefore, lower concentrations of these peptides were observed. Beside the downregulated compounds discussed above, some metabolites were observed at high concentrations, such as nucleotide, base, and proline. Proline has been reported to play important roles in protecting protein denaturation and stabilising protein synthesis, the upregulation of which commonly occurs in cells under heavy metal stresses [43,44]. Because energy generation and central carbon metabolism were severely affected and the upregulation of exonuclease in P.
aeruginosa PAO1 (PA4937) and P. putida KT2440 (PP_0034, PP_0353) was observed (less than 2 fold), the presence of high concentrations of nucleotide and base may be produced by the recycling of denatured nucleic acid.

Conclusions
In this study, we compared the zinc-induced cellular responses of three prominent   [49]. DESeq generated all pairwise comparisons of treatments and associated adjusted P-values (P adj ) controlling for the false discovery rate.
Genes with average absolute fold changes higher than 2.0 and P adj values less than 0.01 were classified as DEGs. The functional clusters of orthologous groups classification of all DEGs was performed automatically using the Integrated Microbial Genomes portal [50].
Validation of RNA-seq results via RT-qPCR To verify the validity of the RNA-seq data, RT-qPCR was performed with the same RNA samples used for RNA-seq analysis. Six genes were selected from each Pseudomonas strain. Therefore, a total of 18 genes were used for validation. The gene-specific primer pairs were designed using Beacon designer (Additonal file 10: and RT-qPCR [25]. Given that all data followed a normal distribution, Pearson's coefficients were calculated to illustrate the correlation level between two data groups.

Metabolite extraction and UPLC-MS/MS analysis
The bacterial samples used for the metabolome analysis were prepared as described by Chavez-Dozal et al. [52]. Zinc-treated Pseudomonas cells were collected from each biological replicate via centrifugation at10,000 × g for 5 min at 4℃. After washing two times with ice-cold phosphate-buffered saline (pH 7.4), the cell pellets were frozen with liquid nitrogen. To extract the metabolites, 50 mg of the frozen samples were dispersed in 1 mL of liquid containing methanol/acetonitrile/water (2:2:1), sonicated for 10 min at 0 ℃, and then centrifuged at 14,000 × g for 20 min at 4 ℃. The resulting supernatants were lyophilised and reconstituted in 100 µL of acetonitrile/water (1:1) before UPLC-MS/MS analysis.
An Agilent 1290 HPLC system (Agilent, Santa Clara, CA, USA) equipped with an HILIC analytical column (amide 1.7 µm, 2.1 mm × 100 mm, Waters, Milford, MA, USA) was used to analyse the extract of each sample. The mobil phase was water (containing 25 mmol L − 1 CH 3 COONH 4 and 25 mmol L − 1 NH 3 .H 2 O, Solvent A), and acetonitrile (Solvent B). The dilution programme was operated for 0.5 min 95:5 (B:A) then switched to 40:60 for 8.5 min before returning to 95:5 for 9 min further. The mass spectra were recorded by an Agilent 6540 Q-TOF mass spectrometer operating in the positive and negative ion modes using an Agilent electrospray ionization source. After each spectrum was normalised to the total spectral intensity, the profile data were imported into SIMCA-P + software (Umetrics AB, Umeå, Sweden) and pareto-scaled to perform orthogonal partial least squares-discriminate analysis (OPLS-DA). The significance of the metabolites was ranked using the variable importance in projection score from the OPLS-DA model. For the univariate analysis of a specific metabolite, statistical significance was determined using Student's two-sample t-test. The metabolites with average absolute fold changes higher than 2.0 and P adj values less than 0.05 were defined as statistically different.

Declarations
Ethic approval and consent to participate Not applicable

Consent for publication
Not applicable Availability of data and materials The RNA-seq datasets generated during the current study are available in the NCBI short read archive database under the bioproject accession number PRJNA606809.

Competing interests
The authors declare that they have no competing interests. Authors' contributions PL and LM conceived the experiment; JC and WL performed the experiment; PL analyzed the transcriptomic and metabolic data obtained in this study; PL wrote the manuscript.