Method for analyzing cadmium stress response of sorghum

Abstract

The invention relates to a sorghum cadmium stress response analysis method, and belongs to the technical field of crop stress response analysis, and the method comprises the following steps: S1, sorghum cadmium stress treatment: selecting uniformly growing sorghum seedlings, growing in a hydroponic culture solution for 7 days, and then treating the seedlings in a CdCl < 2 > solution for 3 days; after treatment, the roots of the sorghum seedlings are cleaned with sterile distilled water, and the sorghum seedlings are immediately frozen in liquid nitrogen; s2, carrying out transcriptome biological information analysis and metabonomics analysis on the processed sorghum roots; and S3, processing the collected data, identifying a differential expression gene of the sorghum root in response to Cd stress, and screening out differential metabolites. According to the analysis method, the key genes and the key metabolites of the sorghum responding to the cadmium stress are excavated.

Description

Technical Field

[0001] The invention belongs to the technical field of crop stress response analysis, and particularly relates to an analysis method for sorghum cadmium stress response. Background Art

[0002] In recent decades, environmental pollution has become increasingly serious, and cadmium (Cd) pollution has seriously threatened plant growth and human health. Cadmium is highly toxic and can severely inhibit the growth and development of plant roots, stems and other organs, leading to a decrease in plant yield and quality.

[0003] Sorghum is a C4 plant that is widely cultivated as a food, feed, and energy crop. Sorghum can be used to remediate Cd-contaminated soils due to its high biomass yield, high tolerance to abiotic stresses, and relatively high Cd accumulation observed in some genotypes. Contaminated plant tissues can be harvested and used to produce ethanol gas, thus preventing Cd contamination from re-entering the food chain. This means that sorghum can serve as a link between phytoremediation and bioethanol production, achieving a win-win effect. Therefore, it is crucial to understand the Cd toxicity of this species and its response to Cd. Summary of the invention

[0004] The purpose of the present invention is to provide an analysis method for cadmium stress response of sorghum and to explore key genes and key metabolites of cadmium stress response of sorghum.

[0005] The technical solution adopted by the present invention to achieve its purpose is as follows:

[0006] A method for analyzing sorghum cadmium stress response comprises the following steps:

[0007] S1. Cadmium stress treatment of sorghum: uniformly grown sorghum seedlings were selected and grown in hydroponic culture medium for 7 days, and then the seedlings were placed in CdCl2 solution for 3 days; after treatment, the roots of the treated sorghum seedlings were washed with sterile distilled water, dried with hand towels, and immediately frozen in liquid nitrogen;

[0008] S2. Perform transcriptome bioinformatics analysis and metabolomics analysis on the treated sorghum roots;

[0009] S3. Process the collected data, identify the differentially expressed genes in sorghum roots in response to Cd stress, and screen out differential metabolites.

[0010] Furthermore, the concentration of the CdCl2 solution in step S1 is 0 μM and 10 μM.

[0011] Furthermore, in step S2, transcriptome bioinformatics analysis includes raw sequencing data, data quality control, comparison with reference genome, quantitative analysis of gene expression, and differential gene screening.

[0012] Furthermore, differential gene screening included GO, KEGG annotation and enrichment analysis, GSEA enrichment analysis, PPI analysis and transcription factor analysis.

[0013] Furthermore, 2683 differentially expressed genes in sorghum roots in response to Cd stress were identified, of which 1812 were up-regulated and 871 were down-regulated.

[0014] Further, in step S2, a UPLC-MS / MS system is used to perform metabolomics analysis, and the metabolomics analysis includes orthogonal partial least squares discriminant analysis and differential metabolite screening. The criteria for differential metabolite screening are FC ≥ 2 and fold change FC ≤ 0.5, and metabolites with variable importance VIP score > 1 in the projection.

[0015] Further, metabolomics analysis includes the following steps:

[0016] A. Sample extraction: The tissue part is the root of sorghum. The sample was taken out from the -80℃ refrigerator and placed in a freeze dryer for vacuum freeze drying. The freeze-dried root was crushed with a grinder at a frequency of 30Hz for 1.5 minutes. 50mg of sample powder was weighed and added to 1200μL of 70% methanol aqueous solution precooled at -20℃, vortexed once for 30min, 30s each time, vortexed 6 times in total, and then left overnight. After the extract was centrifuged at 12000rpm for 3 minutes, the supernatant was aspirated, filtered through a 0.22μm pore size microporous membrane, stored in an injection bottle, and the extract was analyzed using a UPLC-MS / MS system.

[0017] B. Analysis: Analyst 1.6.3 software was used to analyze the metabolite data. Metabolites with fold change FC ≥ 2 and fold change FC ≤ 0.5, and variable importance VIP score in projection > 1 were considered differentially expressed metabolites to screen out differential metabolites.

[0018] Further, the conditions of liquid chromatography are:

[0019] Chromatographic column: Agilent SB-C18 1.8 μm, 2.1 mm*100 mm;

[0020] Mobile phase A: ultrapure water containing 0.1% formic acid; Mobile phase B: acetonitrile containing 0.1% formic acid;

[0021] Instrument column temperature: 40°C; flow rate: 0.35mL / min; injection volume: 2uL;

[0022] Elution gradient: 0.0min, B phase ratio is 5%; 9.0min, B phase ratio increases linearly to 95% and maintains at 95% for 1min, 10.0-11.0min, B phase ratio decreases to 5%; and equilibrates at 5% until 14.0min;

[0023] The mass spectrometry conditions were:

[0024] The temperature of the electrospray ion source was 500°C; the ion spray voltage was 5500 V (positive ion mode) / -4500 V (negative ion mode); the ion source gas I (GSI) was set to 50 psi, the ion source gas II (GSII) was set to 60 psi, the curtain gas (CUR) was set to 25 psi, and the collision induced ionization parameter was set to high; the QQQ scan used the MRM mode, and the collision gas nitrogen was set to medium; through further optimization of the declustering voltage DP and collision energy CE, the DP and CE of each MRM ion pair were completed, and a specific set of MRM ion pairs were monitored in each period according to the metabolites eluted in each period.

[0025] Furthermore, 160 differential metabolites were screened, including 35 down-regulated and 125 up-regulated metabolites.

[0026] Furthermore, the transcriptome bioinformatics analysis was integrated with the metabolome analysis to identify the flavonoid biosynthesis pathway as playing a key role in the response of sorghum roots to Cd stress.

[0027] The beneficial effects of the present invention are:

[0028] Through the analysis method of the present invention, combined with transcriptome and metabolic results, we found that the flavonoid biosynthesis pathway is a key pathway for sorghum seedlings to respond to Cd stress, and multiple DEGs and DEMs in sorghum roots are activated in response to Cd stress. The upregulation of flavonoid biosynthesis pathway genes and the increase of flavonoid metabolites are closely related to sorghum tolerance to Cd stress. These results confirm the importance of the flavonoid metabolic pathway in sorghum root response to Cd stress. The present invention provides new insights into the physiological and molecular mechanisms of sorghum response to Cd stress.

[0029] The present invention identified 2683 differentially expressed genes in sorghum roots in response to Cd stress, of which 1812 were upregulated and 871 were downregulated; and screened out 160 differential metabolites, including 35 downregulated and 125 upregulated metabolites. BRIEF DESCRIPTION OF THE DRAWINGS

[0030] Figure 1 is the sample correlation plot.

[0031] Figure 2 It is the MA map of differentially expressed genes.

[0032] Figure 3 It is the KEGG enrichment bar chart of DEG.

[0033] Figure 4 It is a GO classification bar chart.

[0034] Figure 5 It is a GO enrichment bar chart of differentially expressed genes.

[0035] Figure 6 It is a KOG classification bar chart.

[0036] Figure 7 It is the hierarchical clustering tree diagram of the WGCNA analysis module.

[0037] Figure 8 is the heat map of inter-module correlation.

[0038] Fig. 9 It is the heat map of correlation between samples and modules.

[0039] Fig.10 It is the module gene clustering heat map.

[0040] Fig.11 It is a diagram of the module gene expression pattern.

[0041] Fig.12 is a bar graph of differential metabolites.

[0042] Fig.13 is a volcano plot of differential metabolites.

[0043] Fig.14 It is the KEGG enrichment map of differential metabolites.

[0044] Fig.15 is a plot of differential abundance scores.

[0045] Fig.16 Heat map of flavonoid biosynthesis metabolite expression.

[0046] Fig.17 is the histogram of KEGG enrichment analysis of sorghum root DEGs and DEMs.

[0047] Fig.18 It is the KGML network diagram of DEG and DEM.

[0048] Fig.19 The DEGs and DEMs involved in flavonoid biosynthesis in sorghum roots under Cd stress were shown in parallel. DETAILED DESCRIPTION

[0049] Example 1

[0050] A method for analyzing sorghum cadmium stress response comprises the following steps:

[0051] 1. Sorghum cadmium stress treatment: Sorghum material L407B (provided by the Sorghum Research Laboratory of the Millet Research Institute of Hebei Academy of Agricultural and Forestry Sciences) was used to select healthy sorghum seedlings with uniform growth, grow them in hydroponic culture medium for 7 days, and then place the seedlings in CdCl2 (0 and 10 μM) solution for treatment for 3 days; after treatment, the roots of the treated sorghum seedlings were washed with sterile distilled water, dried with hand towels, and immediately frozen in liquid nitrogen and stored at -80°C for later use;

[0052] 2. Transcriptome bioinformatics analysis and metabolomics analysis of the treated sorghum roots:

[0053] 2.1 Transcriptome bioinformatics analysis:

[0054] The experimental process of transcriptome sequencing includes RNA extraction, RNA detection, library construction and sequencing.

[0055] A. Two ways to obtain mRNA: ① Enrich mRNA with polyA tail using magnetic beads with Oligo (dT); ② Extract ribosomal RNA from total RNA to obtain mRNA; then add fragmentation buffer to break mRNA into short fragments, and use this as a template to synthesize the first-strand cDNA with six-base random primers, then add buffer, dNTPs (dTTP, dATP, dGTP and dCTP) and DNA polymerase I to synthesize the second-strand cDNA, use AMPure XP beads to purify double-stranded cDNA, then perform end repair, add A tail and connect sequencing adapters, screen out cDNA fragments of appropriate length and PCR amplify to obtain the final library. Use Agilent2100 to detect the insert size of the library, use Q-PCR to quantify the effective concentration of the library (>2nM), complete the library picking, and use Illumina platform for sequencing.

[0056] B. Quality control of raw data: filter the raw data to remove reads with adapters; remove paired reads with N content exceeding 10% of the base number of the read; remove paired reads with low quality (Q≤20) base number exceeding 50% of the base number of the read;

[0057] C. Alignment of transcriptome data with reference genome

[0058] HISAT2 software is used to align the original data after quality control with the reference genome to obtain the position information on the reference genome or gene, as well as the sequence feature information unique to the sequencing sample. Based on the position information of the aligned genome, AtringTie is used to assemble the reads into transcripts, and the spliced ​​transcripts are then compared with the annotation information of the genome using GffCompare to discover new transcripts or new genes.

[0059] D. Gene expression quantification: According to the alignment results and the location information of the gene on the reference genome, the number of reads of each gene was counted, and FPKM was used as an indicator to measure the transcript or gene expression level. The calculation formula is: FPKM = number of fragments aligned to a certain transcriptome (number of paired-end reads) / (total number of fragments aligned to the transcriptome × transcript length); the total number of fragments aligned to the transcriptome is in units of 10^6, and the transcript length is in units of 10^3 bases; the FPKM expression of the gene is as follows:

[0060]

[0061] Note: ID represents the gene number, and the second column and subsequent columns represent the FPKM expression levels of all genes in each sample.

[0062] E. Sample correlation analysis: Set up biological replicate experiments and use the Pearson correlation coefficient r as an evaluation indicator of biological replicate correlation. The closer the absolute value of r is to 1, the stronger the correlation between the two replicate samples. Figure 1 As shown;

[0063] F. Differential gene screening: For samples with biological replication, DESeq2 was used to perform differential expression analysis of the samples to obtain the differentially expressed gene set between the two biological conditions. After differential analysis, the Benjamini-Hochberg method was used to perform multiple hypothesis testing correction on the hypothesis test probability to obtain the differential error rate; the conditions for differential gene screening were: |log2Fold Change|≥1, and FDR<0.05; after completing the analysis of differentially expressed genes using DESeq2 / edgeR, the total number of genes in each group, the number of up-regulated genes, and the number of down-regulated genes were counted. The results are as follows:

[0064]

[0065] The differentially expressed genes are listed below:

[0066]

[0067] Note: ID represents gene number; the second column and subsequent columns represent sample expression information

[0068] Constructing a differential gene MA map can intuitively display the overall distribution of gene expression levels and differential folds, see Figure 2 The vertical axis represents the log2FC value, which represents the logarithm of the difference fold; the horizontal axis represents the average gene expression level in the two samples; the red dots represent up-regulated gene expression, the green dots represent down-regulated gene expression, and the blue dots represent no significant difference in gene expression.

[0069] G. KEGG enrichment analysis of differentially expressed genes: Pathway significant enrichment analysis uses pathways in the KEGG database as units and applies hypergeometric test to find pathways that are significantly enriched in differentially expressed genes compared with the entire genome background; the calculation formula for hypergeometric distribution is as follows:

[0070] Where N represents the number of genes with KEGG annotations in all genes, n represents the number of differentially expressed genes in N, M represents the number of genes in a KEGG pathway in N, and m represents the number of differentially expressed genes in a KEGG pathway in M. The KEGG enrichment results are as follows:

[0071] Note: KEGG_pathway: KEGG pathway name; ko_id: KEGG pathway ID; Cluster_frequency: the ratio of the number of differentially expressed genes annotated to the pathway to the number of differentially expressed genes with annotations

[0072] KEGG analysis assigned these DEGs to 120 KEGG pathways. Figure 3 The top 50 are shown. The horizontal axis represents the number of differentially expressed genes annotated to the pathway, and the vertical axis represents the name of the KEGG pathway. The numerical value in the figure represents the number of differentially expressed genes annotated to the pathway. The specific numerical value in the brackets is the ratio of the number of differentially expressed genes annotated to the pathway to the number of background genes annotated to the pathway. The label on the far right represents the category to which the KEGG pathway belongs. Among them, "biosynthesis of secondary metabolites" (ko01110), "glutathione (GSH) metabolism" (ko00480), "biosynthesis of flavonoids" (ko00941), "ATP binding cassette (ABC) transporter" (ko02010), "biosynthesis of phenylpropanoids" (ko00940), "metabolic pathway" (ko01100) and "biosynthesis of flavonoids and flavonols" (ko00944) are significantly enriched.

[0073] H. GO analysis of differentially expressed genes: Taking GO-Term in GO database as the unit, the hypergeometric test was applied to find the GO-Trem significantly enriched in the differentially expressed genes compared with the whole genome background; the enriched terms were made into topGO directed acyclic graph to display the GO nodes (Terms) enriched in differentially expressed genes and their hierarchical relationships; GO analysis showed that DEGs can be mainly divided into three ontologies: biological process (BP), molecular function (MF) and cellular component (CC), such as Figure 4 As shown in the figure, these DEGs can be further divided into 32 functional GO-terms, with the horizontal axis representing the secondary GO terms and the vertical axis representing the number of differentially expressed genes in the GO terms. In the BP group, the most abundant terms were "cellular process", "metabolic process" and "response to stimulus". In the MF category, "binding", "catalytic activity", "transcriptional regulatory activity" and "transport activity" were the terms with the highest expression levels; in the CC process, the most abundant term was "cellular anatomical entity"; the results of GO enrichment analysis are as follows:

[0074]

[0075] Note: Ontology: GO ontology type; ID: GO entry ID; Description: Functional description corresponding to the GO entry; GeneRatio: The ratio of the number of differentially expressed genes annotated to the GO entry to the total number of genes

[0076] Differential gene GO enrichment bar chart ( Figure 5 ), the horizontal axis represents the number of differentially expressed genes annotated to the entry, and the vertical axis represents the name of the GO entry; the numbers in the figure represent the number of differentially expressed genes annotated to the entry, the numbers in brackets represent the specific values ​​of the ratio of the differentially expressed genes annotated to the entry to the background genes annotated to the pathway, and the label on the far right represents the category to which the GO entry belongs.

[0077] I. KOG analysis of differentially expressed genes: Use diamond software to align cDNA sequences to the KOG database, then extract the annotations of the KOG database. According to the annotation results, count the number of differentially expressed genes contained in each KOG functional classification, and use the classification statistics of KOG annotations to draw a bar chart, such as Figure 6 As shown, the horizontal axis represents the functional classification (Code) of KOG ID, and the vertical axis represents the number of differentially expressed genes. Different classifications are represented by different colors. The legend is the Code plus its functional description information.

[0078] J. WGCNA analysis: The WGCNA algorithm is based on high-throughput gene messenger mRNA expression data. It first assumes that the gene network obeys a scale-free distribution and defines the gene co-expression correlation matrix and the adjacency function formed by the gene network. Then, the dissimilarity coefficients of different nodes are calculated and a hierarchical clustering tree is constructed based on them. Different branches of the clustering tree represent different gene modules. The degree of gene co-expression within a module is high, while the degree of gene co-expression in different modules is low.

[0079] Data filtering: First, filter the input fpkm expression level file, and use the varFilter function of the genefilter package in the R language to remove genes with low expression levels in all samples and genes with stable expression levels in all samples to improve the accuracy of network construction. The gene list after filtering is as follows:

[0080]

[0081] Soft threshold selection, set the threshold of the square of the correlation coefficient to 0.85. The higher the square of the correlation coefficient, the closer the network is to the distribution without network scale;

[0082] Model hierarchical clustering: WGCNA analysis constructs a clustering tree based on the correlation between gene expression levels and divides the modules, such as Figure 7 ,A color in the figure corresponds to each gene on the cluster tree belonging to the same module.,If some genes always have similar expression changes in a physiological process or different tissues,,then these genes may be functionally related and can be defined as a module.,For the tree diagram in the upper part, the vertical distance represents the distance between two nodes (genes), and the horizontal distance is meaningless;

[0083] Inter-module correlation heat map: Figure 8 , divided into two parts. The upper part clusters the modules according to the module feature values, and the ordinate represents the degree of difference of the nodes; each row and column in the lower part of the graph represents a module. The darker the color of the block (the redder it is), the stronger the correlation;

[0084] Correlation heat map between samples and modules: Fig. 9 ,If the correlation between a module and the sample is significantly higher than that of other modules, it means that this module has the strongest correlation with the sample;

[0085] Module gene clustering heat map: Fig.10 ,Each dendrogram branch in the figure represents a gene, and the darker the color of each point (white→yellow→red), the stronger the correlation between the two genes corresponding to the row and column;

[0086] Module gene expression pattern: e.g. Fig.11, WGCNA can obtain the module gene expression pattern diagram. The upper part is the clustering heat map of the genes in the module, red is high expression, green is low expression, and the lower part is the expression pattern of the module feature value in different samples. This figure can intuitively show the expression trend of the module genes in different samples;

[0087] Gene list of each module: Connectivity k value, expression information and annotations of seven major databases are added to the gene list of each module obtained by WGCNA analysis; genes with high connectivity k value ranking can be considered as core genes; the gene list of each module is as follows:

[0088]

[0089] geneID represents the gene number; moduleaColors represents the module to which it belongs; kTotal represents the total connectivity of the gene; kWithin represents the connectivity of the gene within the module; kOut represents the connectivity of the gene outside the module, kTotal minus kWithin; kDiff: the difference between kWithin and kOut.

[0090] The relationship between network nodes of each module: The relationship between network nodes of each module is as follows:

[0091]

[0092] fromNode: network node gene 1; toNode: network node gene 2; weight: edge weight of the adjacency matrix, representing the connection strength between two nodes (genes).

[0093] In summary, samples from the 0 and 10 μM Cd treatment groups were selected as research materials for transcriptome and metabolome analysis; a total of 6 libraries (control and Cd-treated plants, each with 3 biological replicates) were generated, and 407,987,342 raw data were sequenced. After removing low-quality reads, a total of 397,496,580 clean reads remained. The Q20, Q30, and GC content values ​​of the clean reads were 97.93-98.32%, 94.14-95.02%, and 54.18-55.42%, respectively, indicating that the transcriptome sequencing results were reliable and could be used for further DEG analysis.

[0094] The percentage of reads that could be mapped to the sorghum reference genome was 91.86-93.62%, of which 89.02-90.76% were uniquely mapped, indicating that the selected reference genome was appropriate. Correlation analysis showed that the treatment groups had appropriate intra-group biological repeatability. Based on the standard cutoff value (|fold change|>2 and corrected p value<0.05), we identified differentially expressed genes (DEGs) in sorghum roots in response to Cd stress. A total of 2683 DEGs (1812 upregulated and 871 downregulated) were identified in the control and Cd treatment.

[0095] 2.2 Metabolomics analysis;

[0096] A. Sample extraction: The tissue part is the root of sorghum. The sample was taken out from the -80℃ refrigerator and placed in a freeze dryer for vacuum freeze drying. The freeze-dried root was crushed with a grinder at a frequency of 30Hz for 1.5 minutes. 50mg of sample powder was weighed and added to 1200μL of 70% methanol aqueous solution precooled at -20℃, vortexed once for 30min, 30s each time, vortexed 6 times in total, and then left overnight. After the extract was centrifuged at 12000rpm for 3 minutes, the supernatant was aspirated, filtered through a 0.22μm pore size microporous membrane, stored in an injection bottle, and the extract was analyzed using a UPLC-MS / MS system.

[0097] The conditions for liquid chromatography were:

[0098] Chromatographic column: Agilent SB-C18 1.8 μm, 2.1 mm*100 mm;

[0099] Mobile phase A: ultrapure water containing 0.1% formic acid; Mobile phase B: acetonitrile containing 0.1% formic acid;

[0100] Instrument column temperature: 40°C; flow rate: 0.35mL / min; injection volume: 2uL;

[0101] Elution gradient: 0.0min, B phase ratio is 5%; 9.0min, B phase ratio increases linearly to 95% and maintains at 95% for 1min, 10.0-11.0min, B phase ratio decreases to 5%; and equilibrates at 5% until 14.0min;

[0102] The mass spectrometry conditions were:

[0103] The temperature of the electrospray ion source was 500°C; the ion spray voltage was 5500 V (positive ion mode) / -4500 V (negative ion mode); the ion source gas I (GSI) was set to 50 psi, the ion source gas II (GSII) was set to 60 psi, the curtain gas (CUR) was set to 25 psi, and the collision induced ionization parameter was set to high; the QQQ scan used the MRM mode, and the collision gas nitrogen was set to medium; the DP and CE of each MRM ion pair were completed by further optimization of the declustering voltage DP and collision energy CE, and a specific set of MRM ion pairs were monitored in each period according to the metabolites eluted in each period;

[0104] B. Analysis: Analyst 1.6.3 software was used to analyze the metabolite data. Metabolites with fold change FC ≥ 2 and fold change FC ≤ 0.5, and variable importance VIP score > 1 in the projection were considered as differentially expressed metabolites. Differential metabolites were screened out. The results of some differential metabolite screening are as follows:

[0105]

[0106] Bar graph of differential metabolites Fig.12 As shown in the figure, the horizontal axis is the log2FC of the differential metabolites, that is, the logarithm of the difference multiple of the differential metabolites with base 2, and the vertical axis is the differential metabolites; red represents up-regulated metabolite content, and green represents down-regulated metabolite content. Fig.13 As shown in the figure, each point represents a metabolite, where green points represent down-regulated differential metabolites, red points represent up-regulated differential metabolites, and gray points represent metabolites that are detected but not significantly different; the horizontal axis represents the logarithm of the difference in the relative content of a metabolite in two groups of samples (log2FC). The larger the absolute value of the horizontal axis, the greater the difference in the relative content of the substance in the two groups of samples; under the VIP+FC screening condition, the vertical axis represents the VIP value, and the larger the vertical axis value, the more significant the difference, and the more reliable the differential metabolites screened; under the VIP+FC+P-value screening condition, the vertical axis represents the difference significance level (-log 10 P-value), the size of the dot represents the VIP value.

[0107] KEGG enrichment analysis of differential metabolites: metabolites were annotated using the KEGG database, and KEGG pathway enrichment analysis was performed based on the differential metabolite results, where RichFactor is the ratio of the number of differential metabolites in the corresponding pathway to the total number of metabolites annotated to the pathway. The larger the ratio, the greater the degree of enrichment. P-value is the hypergeometric test p-value, and the hypergeometric distribution calculation formula is: Among them, N represents the number of metabolites with KEGG annotations in all metabolites, n represents the number of differential metabolites in N, M represents the number of metabolites of a KEGG pathway in N, and m represents the number of differential metabolites of a KEGG pathway in M. The closer the P-value is to 0, the more significant the enrichment is. The top 20 pathways ranked by P-value are selected and displayed from small to large, such as Fig.14 The size of the dots in the figure represents the number of significantly differentially enriched metabolites in the corresponding pathways; the horizontal axis represents the Rich Factor corresponding to each pathway, the vertical axis represents the pathway name (sorted by P-value), and the color of the dots reflects the P-value. The redder the better.

[0108] Analysis of overall changes in KEGG metabolic pathways: The differential abundance score is a pathway-based metabolic change analysis method. The differential abundance score can capture the overall changes of all metabolites in a certain pathway. The calculation formula is:

[0109] Differential abundance score = (number of upregulated differential metabolites in the pathway - number of downregulated differential metabolites in the pathway) / number of all metabolites annotated to the pathway

[0110] Select the top 20 pathways ranked by P-value and display them from small to large. Fig.15 . The ordinate indicates the name of the differential pathway (sorted by P-value), and the abscissa indicates the differential abundance score (DA Score). DA Score reflects the overall changes of all metabolites in the metabolic pathway. A score of 1 indicates that the expression trend of all identified metabolites in the pathway is upregulated, and a score of -1 indicates that the expression trend of all identified metabolites in the pathway is downregulated; the length of the line segment indicates the absolute value of the DA Score, and the size of the dots at the end of the line segment indicates the number of differential metabolites in the pathway. The dots are distributed on the left side of the central axis and the longer the line segment, the more the overall expression of the pathway tends to be downregulated. The dots are distributed on the right side of the central axis and the longer the line segment, the more the overall expression of the pathway tends to be upregulated. The larger the dots, the more metabolites there are; the color of the line segment and the dots reflects the size of the P-value. The closer to red, the smaller the P-value, and the closer to purple, the larger the P-value.

[0111] KEGG pathway enrichment analysis showed that these DEMs could be divided into 47 categories, among which the most important ones were metabolic pathways, biosynthesis of secondary metabolites, followed by flavonoid biosynthesis, ABC transporters and others. The results of this analysis suggested that these metabolites may play an important role in the response of sorghum roots to Cd stress.

[0112] We then constructed a heat map of the expression of flavonoid biosynthesis metabolites. Fig.16, to better visualize the expression patterns of related metabolites in the control and Cd-treated control groups. Hesperidin-7-O-glucoside, butyric acid; 7,3',4'-trihydroxyflavanone, naringenin chalcone; 2',4,4',6'-tetrahydroxychalcone, 5,4'-dihydroxy-7-methoxyflavanone (sakura), taxifolin (dihydroquercetin), phloretin, eriodictyol (5,7,3',4'-tetrahydroxyflavanone), luteolin (5,7,3',4'-tetrahydroxyflavone), 3,5,7-trihydroxyflavone (pine leaf yellow), galangin (3,5,7-trihydroxyflavone), naringenin (5,7,4'-trihydroxyflavone) and apigenin; 4',5,7-trihydroxyflavone, were all found to be involved in the biosynthesis of flavonoids ( Fig.16 ), these results are basically consistent with the transcription results.

[0113] 3. Integrated analysis of transcriptome and metabolome analysis:

[0114] To understand the mechanism of sorghum response to Cd stress, the relationship between DEMs and DEGs was analyzed. The results of KEGG pathway enrichment analysis showed that a large number of differential metabolites and differential genes were enriched in the flavonoid biosynthesis pathway, ABC transporter, flavonoid and flavonol biosynthesis pathway and other pathways ( Fig.17 ).like Fig.18 As shown (squares represent genes or gene products, where red indicates up-regulation and green indicates down-regulation; diamonds represent pathways), co-expression network analysis of DEGs and DEMs in the control and Cd-treated groups showed that these terms were mainly enriched in flavonoid biosynthesis-sorghum (sbi00941), isoflavonoid biosynthesis (map00943), and phenylpropanoid biosynthesis-sorghum (sbi00940 and sbi00260). The results suggest that the flavonoid biosynthesis pathway may play a key role in the response of sorghum roots to Cd stress.

[0115] To further investigate the effects of Cd stress on genes and metabolites related to flavonoid biosynthesis in sorghum roots, DEGs and DEMs related to the flavonoid biosynthesis pathway were analyzed. It was found that key genes in the flavonoid biosynthesis pathway, including naringenin 7-O-methyltransferase (NOMT), flavonoid 3′,5′-hydroxylase (F3'5'H), chalcone synthase (CHS), chalcone-flavonoid isomerase (CHI), and F3M (flavonoid 3'-monooxygenase), were up-regulated in the Cd-treated group. The increase in naringenin was consistent with the up-regulation of genes encoding CHS (chalcone synthase) and CHI (chalcone-flavonoid isomerase). In addition, metabolites related to flavonoid biosynthesis, including hesperidin-7-O-glucoside, accumulated significantly in sorghum under Cd treatment ( Fig.19, triangles and rectangles represent metabolites and genes regulated under Cd treatment, respectively, red indicates up-regulation; blue indicates down-regulation; the types of metabolites are distinguished by long rectangular boxes of different colors, CHS: chalcone synthase; CHI: chalcone-flavonoid isomerase; FLS: flavonol synthase; F3M: flavonoid 3'-monooxygenase; F3'5'H: flavonoid 3',5'-hydroxylase; NOMT: naringenin 7-O-methyltransferase; FLAS: flavonol synthase II).

[0116] The comprehensive analysis results showed that DEGs and DEMs related to the flavonoid biosynthesis pathway significantly responded to Cd stress in sorghum roots.

Claims

1. A method for analyzing cadmium stress response in sorghum, characterized in that: The following steps are involved: S1. Sorghum cadmium stress treatment: uniformly grown sorghum seedlings were selected and grown in hydroponic culture medium for 7 days, and then the seedlings were placed in CdCl2 solution for 3 days; after treatment, the roots of the treated sorghum seedlings were washed with sterile distilled water and immediately frozen in liquid nitrogen; S2. Perform transcriptome bioinformatics analysis and metabolomics analysis on the treated sorghum roots; S3. Process the collected data, identify the differentially expressed genes in sorghum roots in response to Cd stress, and screen out differential metabolites.

2. The method for analyzing cadmium stress response of sorghum according to claim 1, characterized in that: The concentration of the CdCl2 solution in step S1 is 0 μM and 10 μM.

3. The method for analyzing cadmium stress response of sorghum according to claim 1, characterized in that: In step S2, transcriptome bioinformatics analysis includes raw sequencing data, data quality control, comparison with reference genome, quantitative analysis of gene expression, and screening of differentially expressed genes.

4. The method for analyzing cadmium stress response of sorghum according to claim 3, characterized in that: Differential gene screening included GO enrichment analysis, KEGG enrichment analysis, and GSEA enrichment analysis.

5. The method for analyzing cadmium stress response of sorghum according to claim 1, characterized in that: A total of 2683 differentially expressed genes in sorghum roots in response to Cd stress were identified, of which 1812 were up-regulated and 871 were down-regulated.

6. The method for analyzing cadmium stress response of sorghum according to claim 1, characterized in that: In step S2, a UPLC-MS / MS system is used for metabolomics analysis, which includes orthogonal partial least squares discriminant analysis and differential metabolite screening. The criteria for differential metabolite screening are FC ≥ 2 and fold change FC ≤ 0.5, as well as metabolites with variable importance VIP score > 1 in the projection.

7. The method for analyzing sorghum cadmium stress response according to claim 1 or 6, characterized in that: Metabolomics analysis includes the following steps: A. Sample extraction: The tissue part is the root of sorghum. The sample was taken out from the -80℃ refrigerator and placed in a freeze dryer for vacuum freeze drying. The freeze-dried root was crushed with a grinder at a frequency of 30Hz for 1.5 minutes. 50mg of sample powder was weighed and added to 1200μL of 70% methanol aqueous solution precooled at -20℃, vortexed once for 30min, 30s each time, vortexed 6 times in total, and then left overnight. After the extract was centrifuged at 12000rpm for 3 minutes, the supernatant was aspirated, filtered through a 0.22μm pore size microporous membrane, stored in an injection bottle, and the extract was analyzed using a UPLC-MS / MS system. B. Analysis: Analyst 1.6.3 software was used to analyze the metabolite data. Metabolites with fold change FC ≥ 2 and fold change FC ≤ 0.5, and variable importance VIP score in projection > 1 were considered differentially expressed metabolites to screen out differential metabolites.

8. The method for analyzing cadmium stress response of sorghum according to claim 7, characterized in that: The conditions for liquid chromatography were: Chromatographic column: Agilent SB-C18 1.8 μm, 2.1 mm*100 mm; Mobile phase A: ultrapure water containing 0.1% formic acid; Mobile phase B: acetonitrile containing 0.1% formic acid; Instrument column temperature: 40°C; flow rate: 0.35mL / min; injection volume: 2uL; Elution gradient: 0.0 min, B phase ratio 5%; At 9.0 min, the proportion of phase B increased linearly to 95% and was maintained at 95% for 1 min. At 10.0-11.0 min, the proportion of phase B decreased to 5% and was balanced at 5% until 14.0 min. The mass spectrometry conditions were: The temperature of the electrospray ion source was 500°C; the ion spray voltage was 5500 V (positive ion mode) / -4500 V (negative ion mode); the ion source gas I (GSI) was set to 50 psi, the ion source gas II (GSII) was set to 60 psi, the curtain gas (CUR) was set to 25 psi, and the collision induced ionization parameter was set to high; the QQQ scan used the MRM mode, and the collision gas nitrogen was set to medium; through further optimization of the declustering voltage DP and collision energy CE, the DP and CE of each MRM ion pair were completed, and a specific set of MRM ion pairs were monitored in each period according to the metabolites eluted in each period.

9. The method for analyzing cadmium stress response of sorghum according to claim 1, characterized in that: A total of 160 differential metabolites were screened, including 35 down-regulated and 125 up-regulated metabolites.

10. The method for analyzing cadmium stress response of sorghum according to claim 1, characterized in that: The integrated analysis of transcriptome bioinformatics and metabolome analysis identified the flavonoid biosynthesis pathway as playing a key role in the response of sorghum roots to Cd stress.

Images