Integrated Transcriptomic and Metabolomic Analyses of the Molecular Mechanisms of Two Highland Barley Genotypes with Differential Pyroxsulam Responses

Background: Highland barley is one of the few crops that can be grown at high elevations, making it a key resource within the Tibet Plateau. Weeds are a signi�cant threat to highland barley production and new herbicides and tolerant barley varieties are needed to control this ever-growing problem. Results: A better understanding of existing herbicide resistance mechanisms is therefore needed. In this study, transcriptomic and metabolomic analyses were used to identify molecular and physiological changes in two highland barley genotypes with differing sensitivities to pyroxsulam. We identi�ed several stress-responsive metabolites, including �avonoids and antioxidants, which accumulated to signi�cantly higher levels in the pyroxsulam-resistant genotype. Additionally, we found key genes in both the �avonoid biosynthesis pathway and the antioxidant system were up-regulated in pyroxsulam-resistant barley. Conclusions: This work signi�cantly expands on the current understanding of the molecular mechanisms underlying differing pyroxsulam tolerance among barley genotypes and provides several new avenues to explore for breeding or engineering tolerant barley.


Introduction
Highland barley (Hordeum vulgare L.) is an economically important crop and is the only cereal that can mature normally in the short growing seasons often found at high altitudes [1].It is the fourth most important cereal crop worldwide and is primarily utilized for animal feed, brewing, distilling and malting [2].However, highland barley remains a major food source in some areas, such as West Asia and North Africa [3].Although it is grown throughout the world, the majority of its cultivation lies in East Asian countries, such as China, which produces 77% of the total highland barley grown worldwide [4].Highland barley is also consumed as a staple crop in regions of the Qinghai-Tibet Plateau, including Sichuan, Gansu, Qinghai and Tibet.Highland barley is consumed in this region due to its health bene ts, which stem from its metabolites, bioactive carbohydrates, polyphenols, vitamins, phenolic, avonoids and βglucans [5].Additionally, research has indicated that highland barley may possess antihyperglycemic, antihyperlipidemic and anticancer activities [6].Despite these bene ts, there are several challenges in growing barley in the Tibet Plateau, including a high amount of invasive weeds.However, the highland barley is grown in the harsh environment of the Tibet Plateau, the yield and productivity are signi cantly impaired by troublesome weeds.There are many weeds and single herbicides in highland barley eld, which is the main factor restricting the development of barley industry.The development of herbicideresistant barley varieties is critical to increasing the production e ciency of this crop system by enabling the use of more e cient weed control systems.
Pyroxsulam is a new triazolopyrimidine sulfonamide acetolactate synthase (ALS)-inhibiting herbicide developed by Dow AgroSciences (Indianapolis, IN), which provides broad-spectrum control of many annual, biannual and perennial weeds.It has the advantage of high e cacy with low doses and a favorable environmental pro le [7].Pyroxsulam acts by inhibiting ALS, an important target for herbicides, which is an essential enzyme that catalyzes the rst step in the synthesis of the branched-chain amino acids valine, leucine and isoleucine [8].The resulting lack of amino acids in the plant inhibits DNA synthesis, subsequently stopping cell division and causing death in susceptible plants [9].Commercial formulations of pyroxsulam with the herbicide safener cloquintocet allow for selective weed control in cereal crops.Pyroxsulam is currently one of the most widely used herbicides in the world, but little is known about its effects on highland barley.A better understanding of the mechanisms underlying variable pyroxsulam tolerance among highland barley genotypes is critical to e ciently utilizing this herbicide in its production.
The success of breeding herbicide-tolerant varieties relies on the availability of tolerant locally adapted germplasm.Recently, the integration of various omics technologies has proven to be indispensable for determining potential tolerance mechanisms.We performed a comparative analysis of two highland barleys with differing tolerance to pyroxsulam and identi ed substantial differences in their metabolomes and transcriptomes.These differences included key genes in pathways associated with amino acid biosynthesis, reactive oxygen tolerance and avonoids.These results showed that the pro led metabolic alterations in highland barley of a pyroxsulam-resistive and a pyroxsulam-sensitive, and have a reprogramming of the amino acids biosynthesis pathway toward the avonoid pathway after treatment with herbicide pyroxsulam.Our results imply that pyroxsulam-tolerant highland barley varieties induce the expression of avonoids and reactive oxygen species (ROS) scavengers enzymes in order to resist pyroxsulam toxicity.These results provide a good starting point for understanding highland barley pyroxsulam tolerance, although the underlying mechanisms behind these different responses are complex and require additional research.Our results provided deeper insights into the molecular mechanisms and candidate genes for further study and breeding the resistance against herbicide pyroxsulam in highland barley.

Plant material
Qing 0160 and Qing 0305 were collected from Xining City, Qinghai Province, China(N.36.434419,E.101.450759).The samples were identi ed by Liling Jiang and deposited at the National duplicate Genbank for crops, Xining, Qinghai province, China.The specimen accession number were ZDM1651 and ZDM 8091.

Sample preparation and RNA isolation
Two highland barley genotypes were used for this study, including one pyroxsulam-sensitive genotype and one insensitive genotype.All highland barley seedlings were grown in a potting medium under a 11 h light (25°C)/13 h dark (20°C) day/night cycle with 75% humidity.The pyroxsulam herbicide was applied during the two-leaf stage, with the typical eld application concentration of 12.5 g/666.7 m 2 .Leaf samples were harvested at 0, 1 and 6 days after treatment with pyroxsulam, then stored at − 80°C for metabolite pro ling and RNA-sequencing.Total RNA was isolated using an RNA extraction kit from (Sangon, Shanghai, China).

Library preparation for transcriptome sequencing
Total RNA was used as input material for the RNA sample preparations.Brie y, mRNA was puri ed from total RNA using poly-T oligo-attached magnetic beads.Fragmentation was carried out using divalent cations under elevated temperature in First Strand Synthesis Reaction Buffer (5X).First-strand cDNA was synthesized using a random hexamer primer and M-MuLV Reverse Transcriptase, followed by RNaseH RNA digestion.Second-strand cDNA synthesis was subsequently performed using DNA Polymerase I and dNTPs.The remaining overhangs were converted into blunt ends via exonuclease/polymerase activities.After adenylation of the 3' ends of the DNA fragments, adaptors with hairpin loop structures were ligated to prepare for hybridization.In order to select cDNA fragments that were 370-420 bp in length, the library fragments were puri ed with the AMPure XP system (Beckman Coulter, Beverly, USA).Then, PCR was performed with Phusion High-Fidelity DNA polymerase, Universal PCR primers and an index (X) Primer.Finally, PCR products were puri ed using the AMPure XP system and library quality was assessed on the Agilent Bioanalyzer 2100 system.

Sequence data mapping and transcriptome analysis
Raw FASTQ formatted reads were rst processed through in-house Perl scripts to generate clean reads, with adaptor sequences trimmed and reads containing Ns or low-quality bases removed.Scripts were then used to calculate Q20, Q30 and GC content.All the downstream analyses were based on clean data with only high-quality reads retained.
Reference genome and gene model annotation les were downloaded from a public database.The reference genome index was built using Hisat2 v2.0.5 and paired-end clean reads were aligned to the reference genome using Hisat2 v2.0.5.

Differential expression analysis
Differential expression analysis of the two conditions (three biological replicates per condition) was performed using the DESeq2 R package (1.20.0).The resulting P-values were adjusted using the Benjamini and Hochberg's approach for controlling the false discovery rate.Genes with an adjusted Pvalue < 0.05 found by DESeq2 were considered differentially expressed.

GO and KEGG enrichment analysis of differentially expressed genes
Gene Ontology (GO) enrichment analysis of differentially expressed genes (DEGs) was implemented via the clusterPro ler R package.GO terms with corrected P-values less than 0.05 were considered signi cantly enriched.KEGG was employed to determine high-level function of DEGs (http://www.genome.jp/kegg/).We used the clusterPro ler R package to test the statistical enrichment of DEGs in KEGG pathways.

Metabolite extraction
Each sample consisted of 100 mg of tissue, which was ground with liquid nitrogen, followed by resuspension in prechilled 80% methanol and 0.1% formic acid by vortexing.The samples were incubated on ice for 5 min and then centrifuged at 15,000 rpm and 4°C for 5 min.Some of the supernatant was diluted to a nal concentration containing 53% methanol with LC-MS grade water.The samples were subsequently transferred to a fresh Eppendorf tube and were then centrifuged at 15,000 g and 4°C for 10 min.Finally, the supernatant was injected into an LC-MS/MS system for analysis [1].Liquid sample (100 µL) and prechilled methanol (400 µL) were mixed by vortexing [2][3].Samples and prechilled 80% methanol were mixed by vortexing, and then sonicated for 6 min [4][5].

UHPLC-MS/MS analysis
UHPLC-MS/MS analyses were performed using a Vanquish UHPLC system (Thermo Fisher, Germany) coupled with an Orbitrap Q ExactiveTM HF mass spectrometer (Thermo Fisher, Germany) at the Novogene Co., Ltd.(Beijing, China).Samples were injected onto a Hypersil Gold column (100×2.1 mm, 1.9 µm) using a 17 min linear gradient at a ow rate of 0.2 mL/min.The eluents for the positive polarity mode were eluent A (0.1% FA in water) and eluent B (methanol).The eluents for the negative polarity mode were eluent A (5 mM ammonium acetate, pH 9.0) and eluent B (methanol).The solvent gradient was set as follows: 2% B, 1.5 min; 2-100% B, 12.0 min; 100% B, 14.0 min; 100-2% B, 14.1 min; and 2% B, 17 min.Q ExactiveTM HF mass spectrometer was operated in positive/negative polarity mode with spray voltage of 3.2 kV, capillary temperature of 320°C, sheath gas ow rate of 40 arb and aux gas ow rate of 10 arb.

Data processing and metabolite identi cation
The raw data les generated by UHPLC-MS/MS were processed using Compound Discoverer 3.1 software (CD3.1,Thermo Fisher) to perform peak alignment, peak picking and quantitation for each metabolite.The main parameters were set as follows: retention time tolerance, 0.2 minutes; actual mass tolerance, 5 ppm; signal intensity tolerance, 30%; signal/noise ratio, 3; and minimum intensity, 100,000.Next, peak intensities were normalized to the total spectral intensity.The normalized data were used to predict the molecular formula based on additive ions, molecular ion peaks and fragment ions.Peaks were then matched with the mzCloud (https://www.mzcloud.org/),mzVault and MassList databases to obtain accurate qualitative and relative quantitative results.Statistical analyses were performed using the statistical software R (R version R-3.4.3),Python (Python 2.7.6 version) and CentOS (CentOS release 6.6).
When data were not normally distributed, normal transformations were attempted using the area normalization method.

Global transcriptomic difference in highland barley
In order to determine the different molecular mechanisms of two highland barley genotypes in response to pyroxsulam, we compared transcriptomic data at 0, 1 and 6 days after pyroxsulam treatment.Raw and clean reads were obtained from 18 samples using the Illumina sequencing platform (Table S1).We applied fragments per kilobase of transcript per million mapped reads (FPKM) to normalize the expression levels of the genes.The numbers of DEGs in highland barley with more than a 2-fold change in expression (P < 0.05) are shown in Table S2.The numbers of DEGs varied between treatment-by-time (Fig. 1A) and the volcano maps for all DEGs are shown in Fig. 1B-1D.Among the DEGs, the pyroxsulam without treatment resulted in the largest number of DEGs (6432 up-regulated and 6143 down-regulated genes in R0 vs S0) (Fig. 1B).One day after pyroxsulam treatment resulted in the smallest number of DEGs (2013 up-regulated and 2191 down-regulated genes in R1 vs S1) (Fig. 1C).The number of DEGs after 6 days of pyroxsulam exposure was slightly higher than 0 days after treatment (3967 up-regulated and 3869 down-regulated genes in R6 vs S6) (Fig. 1D).These results suggestion that pyroxsulam different treatment time showed different gene expression pro les.

Gene ontology (GO) enrichment analysis of DEGs
DEGs were subjected to GO enrichment analysis to identify enriched biological functions.Genes associated with small molecule metabolic process, carbohydrate metabolic process and organic acid metabolic process showed the greatest differential expression in the R0vsS0 group in the biological process (BP) category, while genes associated with NAD binding and oxidoreductase activity, acting on the aldehyde or oxo group of donors showed the greatest differential expression in the R0vsS0 group in the molecular function (MF) category.Genes associated with transferase activity, transferring glycosyl groups, transmembrane transporter activity and transferase activity, and transferring hexosyl groups showed the greatest differential expression in the R1vsS1 group in the MF category.Genes associated with small molecule metabolic process, cellular amide metabolic process and peptide metabolic process showed the greatest differential expression in the R6vsS6 group in the BP category.Genes associated with non-membrane-bounded organelle, intracellular non-membrane-bounded organelle and ribonucleoprotein complex showed the greatest differential expression in the R6vsS6 group in the cellular component (CC) category.Finally, genes associated with structural molecule activity, coenzyme binding and structural constituent of ribosome showed the greatest differential expression in the R6vsS6 group in the MF category (Fig. 2).These data indicate that there are signi cant differences in highland barley resistant and sensitive varieties in a variety of biological processes, molecular functions and cellular components in response to pyroxsulam, which likely drive their different tolerances to this herbicide.

KEGG Pathway classi cation of DEGs.
To identify biological pathways related to DEGs, KEGG pathway analysis was performed.We found that the DEGs were associated with several enriched signaling pathways, including "Carbon metabolism", "Biosynthesis of amino acids" and "Oxidative phosphorylation", when comparing R0 to S0 (Fig. 3A).In the R1vsS1 group, the DEGs were associated with "Ribosome" and "Photosynthesis -antenna proteins" (Fig. 3B).In the R6vsS6 group, "Ribosome", "Biosynthesis of amino acids" and "Carbon metabolism" were signi cantly enriched (Fig. 3C).These results indicated that metabolism, biosynthesis, modi cation and other pathways may be involved in these two genotypes' differing tolerance to pyroxsulam.These gene expression differences highlight common responses to the pyroxsulam condition but also show that these genotypes differ substantially from each other.

Screening of differential metabolites of highland barley
Metabolomics is closer to phenotypic omics, and can be viewed as an extension of transcriptomics and proteomics that re ects the current physiological state of an organism.To compare the metabolic changes that took place after pyroxsulam treatment in tolerant and susceptible genotypes, we employed metabolomics.In total, 837 and 386 metabolites in the positive and negative ionization modes were detected and identi ed in all samples (Table S1 and S2).Volcano plots can be used to display the overall distribution of different metabolites, with the abscissa representing the differences in the changes of metabolites in different groups (log 2 FoldChange), and the ordinate representing the signi cance level of differences (-log 10 p-value).Each point in the volcano chart represents a different metabolite.The signi cantly up-regulated metabolites are represented by red dots, while the signi cantly down-regulated metabolites are represented by green dots, and the size of the dot represents the VIP value.After constructing a volcano plot, it could be seen that signi cantly more stress-responsive metabolites, including avonoids and antioxidants, were identi ed in the pyroxsulam-resistant genotype compared to the pyroxsulam-sensitive genotype (Fig. 4).

Differential metabolite KEGG enrichment
To understand the metabolic activities of the two highland barleys after pyroxsulam treatment, we mapped the annotated metabolites to KEGG pathways.Enrichment analysis of the KEGG pathways can determine the main biological functions performed by the differential metabolites.According to the above enrichment results, a bubble chart of the top 20 enriched KEGG pathways was constructed (Fig. 5).A number of pathways were found to be differentially expressed, including "carbon metabolism", "Phenylalanine metabolism", "ABC transporters" and "Biosynthesis of amino acids" in the R0vsS0 comparison group (Figs.5A and B)."Metabolic pathways", "Tryptophan metabolism" and "Flavonoid biosynthesis" were found to have the greatest enrichment in the R1vsS1 group (Figs.5C and D).Finally, "Biosynthesis of secondary metabolites" and "Flavonoid biosynthesis" were most signi cantly enriched in the R6vsS6 group (Figs.5E and 5F).These metabolites cover most of central metabolism and re ect the physiological state.Taken together, these results show that avonoid biosynthesis was consistently different in tolerant versus susceptible highland barley genotypes, indicating that it likely plays a central role in the differing tolerances of the two highland barley varieties.

Correlation between differential gene and differential metabolite expression
To understand the regulatory networks of the two highland barley genotypes' responses to pyroxsulam, a correlation analysis was carried out with DEGs obtained by transcriptome analysis and the signi cantly different metabolites obtained by metabolomics analysis to construct an integrated metabolic map (Fig. 6).The map showed that the pathways involved in amino acid metabolism were enriched signi cantly without pyroxsulam treatment (R0 vs. S0) (Fig. 6A).The gene expression and metabolic ux shifted from amino acid metabolism to pyruvate metabolism and tropane, piperidine and pyridine alkaloid biosynthesis soon after pyroxsulam treatment (R1 vs. S1) (Fig. 6B).The genes involved in the phenylpropanoid and glutathione metabolic pathways were up-regulated, leading to the high accumulation of avonoid and antioxidants after 6 days of pyroxsulam treatment (R6 vs. S6) (Fig. 6C).
Overall, these results suggest that several genes encoding key enzymes involved in the biosynthesis process of avonoids and other crucial metabolites were differentially expressed, suggesting that these genes could be valuable targets for improved pyroxsulam tolerance in highland barley.

Discussion
Weeds pose a signi cant threat to the production of highland barley, and new herbicides are required to deal with this ever-evolving problem [10].Advances in biotechnology have led to the discovery of several herbicides that belong to different chemical classes and possess diverse modes of action.These newer herbicides represent the most cost-effective and environmentally sustainable way to control weeds and increase worldwide food production [11].Previously, herbicides could only be applied to a limited number of crops but the development of herbicide-resistant varieties has greatly expanded their spectrum of use [12,13].Highland barley has not been extensively investigated for herbicide tolerance, which limits the herbicides that can be used to control weeds, thus limiting barley productivity.In this study, we utilized time-course metabolic pro ling and transcriptional pro ling of tolerant and sensitive highland barley accessions to better understand the mechanisms underlying their tolerance differences.An integrated analysis of the transcriptome and metabolome revealed signi cantly more stress-responsive metabolites in the pyroxsulam-resistant highland barley compared to pyroxsulam-susceptible.The up-regulation of key genes in the phenylpropanoid pathway led to avonoid accumulation, which may represent a mechanism responsible for their differing tolerances.In addition, several unique metabolites and unidenti ed compounds, such as glutathione, oxyresveratrol, artemtherin and others, could play a protective role during responses to pyroxsulam.Metabolic reprogramming of the phenylpropanoid pathway has been reported to play a key role in response to biotic and abiotic stresses in rice, maize and soybean [14][15][16].Metabolic pro ling during fungal infection has shown a mobilization of carbohydrates, changes in amino acid pools, and the activation of avonoid biosynthetic pathways in soybean [17].Additionally, drought-stressed Arabidopsis has been reported to possess signi cantly different metabolite pro les, including large changes in avonoids [18].
It has previously been reported that avonoids play a major role in protection against stresses in plants.
For example, avonoids accumulate in maize plants growing at high altitudes to prevent damage caused by high UV-B exposure [19].An increase in avonoid secondary metabolites has also been shown to provide protection during pathogen infection [20].Researchers have also reported that avonoids are involved in the resistance to aluminum toxicity in maize [21].Flavonoids have been shown to possess antimicrobial activity in plants, and have potent antifungal activity against several major plant fungal pathogens, including Fusarium oxysporum [22].In the current study, we found that the biosynthesis of avonoids signi cantly increased during pyroxsulam treatment, with pyroxsulam-tolerant barley showing a stronger induction.These results provide a starting point for using metabolomics to assist in the breeding of pyroxsulam-tolerant highland barley varieties.
Herbicides have been demonstrated to result in the generation of ROS in plants, with antioxidant systems often being responsible for the tolerance level of a given species [23].Thioredoxin reductase plays a key role in catalyzing the NADPH-dependent reduction of oxidized thioredoxin to scavenge ROS, and this gene has been shown to be induced during salt stress in hulless barley [24].Gene expression analysis in wheat revealed that most unigenes encoding peroxidase were induced by biotic stress, leading to better resistance to Bipolaris sorokiniana [25].In the current study, transcriptome analysis revealed that the gene expression levels of glutathione S-transferase and glutathione synthetase, which both scavenge ROS, were signi cantly up-regulated in the pyroxsulam-resistant highland barley genotype.In addition, metabolite analysis revealed that the metabolite glutathione is uniquely induced in the tolerant highland barley after treatment with pyroxsulam.These results demonstrate that genes associated with the antioxidant system were up-regulated, which led to accumulation of antioxidants and subsequent scavenging of ROS in the tolerant highland barley variety.Our results shed new light on the molecular mechanisms underlying variable pyroxsulam tolerance in highland barley, and are an important rst step in the development of tolerant crops through gene editing or traditional breeding.The research to promote sustainability and productivity of an important economic crops highland barley and more e ciently develop resistant crops through gene modifcation, molecular breeding, or novel genetic-editing technologies.

Declarations Figures
The    Each point in the volcano plot represents a metabolite, and the signi cantly up-regulated metabolites are marked with red dots, while signi cantly down-regulated metabolites are represented by green dots.The size of the dot represents the VIP value.
Differential metabolite KEGG enrichment.The left frame is the positive ion mode and the right frame is the negative ion mode.The abscissa is the ratio of the number of differential metabolites in the corresponding metabolic pathway to the total number of metabolites identi ed in the pathway.The larger number of DEGs down-and upregulated in the resistant and sensitive highland barley genotypes after pyroxsulam treatment.A. DEG counts of highland barley transcriptomes.B, C and D. Volcano maps of DEGs in the R0vsS0, R1vsS1 and R6vsS6 groups, respectively.

Figure 2 GO
Figure 2