Gene coexpression network analysis combined with metabonomics reveals the resistance responses to powdery mildew in Tibetan hulless barley

Powdery mildew is a fungal disease that represents a ubiquitous threat to crop plants. Transcriptomic and metabolomic analyses were used to identify molecular and physiological changes in Tibetan hulless barley in response to powdery mildew. There were 3418 genes and 405 metabolites differentially expressed between the complete resistance cultivar G7 and the sensitive cultivar Z13. Weighted gene coexpression network analysis was carried out, and the differentially expressed genes were enriched in five and four major network modules in G7 and Z13, respectively. Further analyses showed that phytohormones, photosynthesis, phenylpropanoid biosynthesis, and flavonoid biosynthesis pathways were altered during Qingke-Blumeria graminis (DC.) f.sp. hordei (Bgh) interaction. Comparative analyses showed a correspondence between gene expression and metabolite profiles, and the activated defenses resulted in changes of metabolites involved in plant defense response, such as phytohormones, lipids, flavone and flavonoids, phenolamides, and phenylpropanoids. This study enabled the identification of Bgh responsive genes and provided new insights into the dynamic physiological changes that occur in Qingke during response to powdery mildew. These findings greatly improve our understanding of the mechanisms of induced defense response in Qingke and will provide new clues for the development of resistant Tibetan hulless barley varieties.

Multi-omics techniques including genomics, transcriptomics, proteomics and metabolomics can be used to generate massive and complex data sets 11 . In combination, a systems biology approach has the potential to comprehensively investigate the complex biological processes that occur when plants respond to powdery mildew and define how these components dynamically interact. Therefore, a better understanding of the reconstruction of biochemical networks 12 and complex signaling pathways can be obtained by the integration and analysis of omics data 13 .
Tibetan hulless barley (Hordeum vulgare L. var. nudum), also called "Qingke" in Chinese and "Ne" in Tibetan, with highland adaptation to extreme environmental conditions, is the principal food for Tibetans and an important crop in the Tibetan Plateau 14 . One of the major diseases of Tibetan hulless barley is powdery mildew caused by Bgh, and the cultivation of Bgh-resistant crop cultivars is the most economical method for controlling the damage. A previous study showed that Qingke varieties exhibited different responses to powdery mildew, suggesting that there is extensive genetic variation to explore 15 . Our limited understanding of Tibetan hulless barley genetics seriously hindered the systematical investigation of genes and molecular mechanisms underlying its resistance response to powdery mildew. Fortunately, using a whole-genome shotgun strategy, we built a 3.89-Gb draft genome of Tibetan hulless barley and predicted 36,151 protein-coding genes to better understand its adaptation to various stressful environments on highland and to facilitate crop improvement 14 . Furthermore, a set of genes involved in plant environmental responses and adaptation, such as plant hormone signal transduction and plant-pathogen interaction, were found to be positively selected in Qingke. However, the gene networks and molecular interactions involved in Tibetan hulless barley resistance to Bgh remain largely unknown. To efficiently breed resistant varieties, it is necessary to better understand how this plant responds to powdery mildew and then use this knowledge to improve the ability to cope with it.
An important object of gene expression studies is to reveal the underlying transcriptional networks, genes and pathways that are mediating a physiological process 16 . Gene coexpression network (GCN) analysis is a powerful systems biology approach whereby modules of highly coexpressed or connected genes can be identified and provide a meaningful way to examine the correlations in gene expression generated from complex RNA-seq datasets across developmental stages, treatments and time courses. A GCN is a set of relationships between genes where a node is defined as a gene connected to other genes by edges based on pairwise similarities 16 . A weighted gene coexpression network approach (WGCNA) assigns weights to the edges based on the strength of the correlation and can be used for finding modules of highly correlated genes, for summarizing clusters using the module eigengene or an intramodular hub gene 17 . WGCNA analysis of 10 Arabidopsis ecotypes following cold, heat, high-light, salt and flagellin treatment either separately or in combination identified several stress regulatory modules as important mediators of the regulatory response 18 . One of the important goals of network analysis is to connect gene expression data to the trait or stress response phenotype. Hub nodes have been found to play vital roles in many networks. Highly connected hub genes are expected to play an important role in biology but not always be significantly associated with the trait of interest 19 . It has been suggested that intramodular hub genes that are highly connected within a module are more likely to be biologically significant if that module is associated with the trait 19,20 .
Plants can biosynthesize specialized metabolites to adapt to various stresses such as biotic and abiotic stresses 21 . It has been suggested that indole metabolite, jasmonic acid (JA) 22 , salicylic acid (SA) 23 , callose 3 , phytoalexins, resveratrol and phenolic metabolites 24 , polyamine 25 , phenylpropanoids and flavonoids are all important compounds associated with plant responses to powdery mildew. To understand the underlying mechanism of crops' stress responses and resistances, researchers have focused on the signaling perception, transcriptional regulation, functional protein expression and biosynthesis of specialized products in plants. Multi-omics analysis has become a powerful approach to identify gene-to-metabolite correlations 26 . The integration of modern omics technologies provides a great opportunity for better understanding plant defense mechanisms against powdery mildew at molecular and cellular levels. Previous studies have suggested that plant response to powdery mildew is a dynamic process, and that gene expression and metabolite accumulation are temporally regulated 7 . This information suggests that it is essential to investigate the dynamic changes at transcriptional, proteomic and metabolic levels during the plant response to powdery mildew. Transcriptomic studies can predict changes in gene expression, while metabolomic studies can investigate altered functions triggered by these genes or proteins. Therefore, the combination of transcriptomic and metabolic approaches is essential for a deeper understanding of plant responses to powdery mildew.
Here, to identify the changes that occur in the crop Tibetan hulless barley in responses to powdery mildew at molecular and physiological levels, we applied a transcriptomic coexpression network approach combined with metabolic analysis to associate gene expression with physiological responses. Coupled analysis of transcriptome and metabonomic parameters was performed over a long-time course (168 h) in a complete resistance cultivar Gannongda7 (G7) and a sensitive cultivar Zangqing13 (Z13). There were 3418 genes and 405 metabolites differently regulated during the plants' interactions with Bgh between the two cultivars. WGCNA analysis showed that the differentially expressed genes were enriched in five and four major network modules in G7 and Z13, respectively. Further analyses showed that phytohormones, photosynthesis, phenylpropanoid biosynthesis and flavonoid biosynthesis pathways were altered during Qingke-powdery mildew interaction. Comparative analyses showed a correspondence between gene expression patterns and metabolite profiles, and the activated defenses resulted in changes of metabolites involved in plant defense response, including phytohormones, lipids, flavone and flavonoids, phenolamides, and phenylpropanoids. In addition, characterization of the transcription factors demonstrated that the WRKY transcription factor family exhibited a different regulation mechanism in the two cultivars and may play an important role in the response to powdery mildew. Gene coexpression network analysis combined with metabonomics in this study enabled the identification of powdery mildew responsive genes based on their differential expression profiles and provided new insights into the dynamic physiological changes that occur during the continuous process of Tibetan hulless barley response to powdery mildew. Results RNA sequencing and differential expression analysis. Powdery mildew can be observed as small grayish patches of fluffy fungal growth on the upper surface of the leaves. These spots resemble small cushions of white powder. The fungus only infects the epidermal layer of the leaf, and the tissue on the opposite side of an infected leaf turns pale green to yellow. Leaves remain green and active for some time following infection and then gradually become chlorotic and die off. G7 is a Qingke cultivar with complete resistance to powdery mildew, while Z13 is a sensitive cultivar. At the two-leaf and one needle stage, the well-developed seedlings were selected and inoculated with Bgh. Subsequently, leaves from each individual plant were collected at 0, 6, 36, 72 and 168 h post inoculation (hpi).
The overall number of differentially expressed genes in Qingke G7 increased at 6, 36, and 72 hpi over time but suddenly reduced at 168 hpi (Fig. 1A). The pattern of differential gene expression in Qingke Z13 was strikingly different with the G7 cultivar, especially at 6 and 168 hpi. Taking the time-point of 168 hpi for example, 2618 and 281 genes were significantly up or down-regulated in Z13, respectively, while only 285 and 247 genes were significantly up or down-regulated in G7.
Significantly different responses to Bgh in Qingke between G7 and Z13 cultivars were observed in the heatmap depicting hierarchical clustering of the gene expression data (Fig. 1B,C). In Qingke G7, the expression level of the majority of differentially expressed genes at 168 hpi was much similar with these genes at 0 hpi (Fig. 1B); while in Z13, the directionality of gene differential expression at 168 hpi differed strikingly with that at 0 hpi (Fig. 1C). The similarity between G7 and Z13 in response to Bgh can also be seen in the heatmap that a set of genes are expressed similarly at 36 hpi and 72 hpi.
A set of genes with similar expression patterns are often functionally correlated. The trends of all the differentially expressed genes in G7 and Z13 were analyzed respectively by Short Time-series Expression Miner software (STEM). A total of 20 trends were obtained from the expression data (  1E). Surprisingly, 372 genes and 335 genes were further found in common in profile18 and profile1 in G7 and Z13, respectively. This result indicated that a great number of common genes (accounting for a large proportion) with similar expression patterns were involved in the process of Qingke-powdery mildew interaction in the two different Qingke varieties.
Interestingly, 1199 (31.4%) up-regulated genes were in common between G7 and Z13 in their response to powdery mildew (Fig. 2C), and 526 (42%) down-regulated genes (account for 57.6% and 60.8% of the down-regulated genes in G7 and Z13, respectively) were in common between the two varieties ( Fig. 2D). This result indicated that the two Qingke varieties triggered many of the same sets of genes and might have a similar resistant mechanism during the process of plant-powdery mildew interaction.  WGCNA based on data collected on the 2520 (G7) and 4147 (Z13) differentially expressed genes, respectively. GCNs are composed of genes that have similar profiles and are highly correlated with each other. The weighted interaction network is shown in Fig. 3A,B (Data S1). Nodes (genes) are connected by edges (coexpression relationships). The connection between two nodes was determined by the correlation between the expression levels of the genes those nodes represent across all experiments used in the analysis.

Coexpression Network
This analysis resulted in a network of G7 that grouped 2520 genes into 8 modules (Fig. S1), the most strongly interconnected of which are shown in Fig. 3A. Here, 5 unique modules were identified, with each module depicted with a different color. The module's gene expression profile was represented by its eigengene (Fig. 3C), its most notable component. Notably, MEbrown comprised genes that were highly expressed (upregulated) at 0 hpi and 168 hpi, and MEblue included highly expressed genes at 36 and 72 hpi. MEblack were highly expressed uniquely at 168 hpi. For example, 164 genes involved in the MEblack module were highly expressed specifically at 168 hpi, which might indicate that this set of genes might be responsible for the process of resistance/adaption of G7 to powdery mildew. A total of 795 genes involved in the MEblue module were highly accumulated, specifically at 6, 36, and 72 hpi, which showed that this set of genes in the MEblue module might be responsible for the process of the early resistance/adaption of G7 to powdery mildew.
In Z13, a network that grouped 4147 genes into 4 modules (Fig. 3B) were identified, with each module depicted with a different color. The module's gene expression profile was represented by its eigengene (Fig. 3D). Notably, MEgreen comprised genes that were highly expressed only at 0 hpi, with those of MEblue uniquely expressed at 168 hpi. Interestingly, a total of 1848 genes involved in the MEblue module were significantly upregulated specifically at 168 hpi, which might indicate that this group of genes might play important roles in the responses of Z13 to powdery mildew.  Investigation of major modules correlated with Tibetan hulless barley response to powdery mildew. As seen in the dendrogram on Fig. 3, each module's gene expression profile is represented by its eigengene. In G7, the five resulting eigengenes correlate with the time points due to their gene expression profiles (Fig. 4A). Notably, 3 modules have similar expression patterns in the three time points (6 h, 36 h, and 72 h), including MEyellow, MEblue and MEpink. While MEblack was upregulated at 168 hpi, and MEbrown comprised genes that were upregulated at both 0 hpi and 168 hpi. Therefore, each of these modules identifies a set of genes of G7 in response to powdery mildew at specific time points.
For example, the MEbrown module had the highest correlation with G7, in which 350 genes were highly specifically expressed in G7 at 0 hpi and 168 hpi, which indicated that this set of genes might be responsible for the resistance of G7 to powdery mildew. KEGG pathway enrichment analysis showed that these genes were significantly enriched in phenylalanine metabolism, terpenoid backbone biosynthesis, sesquiterpenoid and triterpenoid biosynthesis, flavonoid biosynthesis, ubiquinone and other terpenoid-quinone biosynthesis, circadian rhythm-plant, zeatin biosynthesis, isoflavonoid biosynthesis and so on (Table S4).
As seen in the dendrogram on Fig. 5A, each of these modules (MEturquoise, MEgreen, MEblue and MEyellow) identifies a set of genes in Z13 at specific time points. Notably, the MEturquoise module comprising genes that were upregulated in Z13 at 6, 36 and 72 hpi, had the highest correlation (cor = 0.41, P = 3.1e-56) with Z13 in response to powdery mildew. MEblue and MEyellow comprised genes that were upregulated in Z13 at 168 hpi, while the MEgreen module comprised genes that were upregulated at the 0 hpi and 168 hpi.
As an example, 679 genes in the MEturquoise module were specifically enriched in photosynthesis, phenylpropanoid biosynthesis, circadian rhythm-plant, monoterpenoid biosynthesis, plant-pathogen interaction, plant hormone signal transduction, starch and sucrose metabolism, and flavonoid biosynthesis (Table S5), which indicated that this group of genes might be responsible for the responses of Z13 to powdery mildew. A set of genes involved in the MEblue and MEyellow modules were specifically upregulated at 168 hpi and enriched in proteasome, protein processing in endoplasmic reticulum, ABC transporters, pyruvate metabolism, citrate cycle (TCA cycle), glutathione metabolism and ribosome-related biological processes, indicating that these groups of genes might be involved in the later process of powdery mildew infection.

Identification of differentially expressed Transcription factors (TFs).
In consideration of the important regulatory function of TFs in response to various stresses, we analyzed TFs-encoding genes by blast against the Plant Transcription Factor Database (PlnTFDB,V3.0) (http://plntfdb.bio.uni-potsdam.de/v3.0/) 28 . As a result, 275 differentially expressed TFs distributed in 39 families were identified (Table S6). These TFs mainly include the following families: MYB-related (24 genes), bZIP (basic region/leucine zipper motif) (23 genes), NAC (NAM, ATAF1-2, and CUC2) (22 genes), WRKY (18 genes), bHLH (basic helix-loop-helix) (17 genes), MYB (myeloblastosis) (17 genes), and Co-like (14 genes). More details of the identified differentially expressed TFs are provided in Table S6. The genes belonging to WRKY families were further analyzed because of their important roles in plant resistant response to pathogenic bacteria. As a result, 10 KRKY genes were differentially expressed in G7 and Z13, most of which were highly induced in G7 at 36 hpi, while in Z13 they were highly induced at 168 hpi. This result reveals the different transcriptomic regulation mechanism in Qingke varieties in response to powdery mildew.
Widely targeted metabolome analysis. To further analyze the dynamic changes at the metabolite level of Tibetan hulless barley caused by gene expression regulation, widely targeted metabolome analysis was carried out. As a result, a total of 568 known and unknown metabolites were detected and quantified in Tibetan hulless barley during the interaction between plants and powdery mildew (Table S7). By mapping the Metware metabolite database (local database) and the general biochemical pathways based on KEGG, the metabolites were divided into several classes, including amino acids and their derivatives, lipids, nucleotides and their derivates, organic acid and its derivatives, flavone and flavonoids, phenolamides, phenylpropanoids, terpenoids, and phytohormones. (Table S7).
Principal component analysis (PCA) was further used to preliminary understand the overall metabolic differences between samples and variation between groups. The results show that there were differences among the Tibetan hulless barley samples at different time points and that there was an obvious variation tendency among the two Qingke groups, especially in Z13-0 hpi (Mcw01, Mcw02, Mcw03 (Fig. S2). Interestingly, metabolic analyses showed that the levels of most amino acids, phenolamides, and lipids including glycerophospholipid and fatty acid were much higher in G7 than those in Z13 at the 0 hpi time-point. However, at the 168 hpi time-point, the levels of amino acids, amino acid derivatives, and phenolamides were much lower in G7 than those in Z13. The flavone and flavonoids were differentially accumulated at 0 hpi and 168 hpi in G7 and in Z13. The heatmap of metabolites detected in G7 and Z13 are shown in Fig. S3.
Jasmonates are critical components in mediating plant stress-induced systemic signals to activate defense-related genes. We found that differential expression of genes induced by Bgh in Qingke were correlated with changes at the physiological level. Our data showed that a set of JA-related genes were induced by Bgh at different time points, leading to the metabolite changes of jasmonates in G7 and Z13. (Fig. 6B).
At 0 hpi (un-inoculated control), the JA and JA-ILE levels were significantly higher in G7 (1.69 ng/mL and 1.47 ng/mL, respectively) than those in Z13 (0.64 ng/mL and 0.39 ng/mL, respectively), while MEJA was almost the same in the two cultivars. After inoculation with Bgh, endogenous JAs accumulated rapidly in G7 and reached the maxima at 6 hpi (2.78 ng/mL). Then, the JA content recovered to almost the same as that of the un-inoculated control at 36 hpi and 72 hpi and later decreased rapidly to 0.72 ng/mL at 168 hpi. However, in Z13, the concentration of JA accumulated to 1.11 ng/mL at 6 hpi and then changed slightly at 36, 72 and 168 hpi. Interestingly, the concentration of JA-ILE decreased gradually during the interaction in G7 with Bgh, from 1.47 ng/mL at 0 hpi to 0.58 ng/mL at 168 hpi, while in Z13, JA-ILE accumulated gradually from 0.39 ng/mL at 0 hpi to 1.49 ng/mL at 168 hpi.

Integrated analyses of the transcriptomic and metabolic datasets.
To gain a deep understanding of the resistant mechanism of G7 and Z13 based on their obvious resistance ability, we further analyzed the gene expression profile combined with the metabolomics analysis at the 0 hpi time point. We found that 412 and 348 genes were lower or higher expressed in Z13 compared with those in G7, and 71 and 59 compounds were differentially accumulated at the metabolic level. The enrichment analysis indicated that these DEGs were significantly enriched in glyoxylate and dicarboxylate metabolism, phenylpropanoid biosynthesis, monoterpenoid biosynthesis, alpha-linolenic acid metabolism, photosynthesis -antenna proteins, flavonoid biosynthesis, and starch and sucrose metabolism, among others (Table S8). The metabolite profiling showed that the differential metabolites mainly included phenylpropanoids and their derivatives, phenolamides, flavonoids, lipids, terpenoids and so on. Further correlation analysis of differential metabolites and DEGs at 0, 6, 36, 72, 168 hpi were carried out. The results showed that a great number of metabolites accumulated differentially during the process of Qingke-Bgh interaction, which were synchronous with the expression of the differentially expressed genes. 443 DEGs were identified to be related to the differential accumulation of 110 metabolites in the process of qingke response to powdery mildew (Table S9). This result indicated that the alteration in gene expression induced by powdery mildew was correlated with the metabolic changes at the physiological level during the process of Qingke-powdery mildew interaction.

Discussion
To promote sustainability and productivity of plants in the face of powdery mildew, we must more efficiently develop resistant crops through gene modification, molecular breeding, or novel genetic-editing technologies. To effectively utilize these strategies to produce the next generation of crops, the multi-omics datasets should be appropriately exploited. The present study reports the first effort to integrate transcriptomic and metabolic techniques for the comparative analyses of the genes and the metabolites involved in responses of Qingke plants to powdery mildew. The results enhance our understanding of the mechanisms underlying the responses of Tibetan hulless barley to powdery mildew.
Previous studies had shown that pathogens cause changes in gene expression that fluctuate over time. To investigate this further, WGCNA analyses were performed in the present study and revealed extensive temporal regulation of transcript levels in Tibetan hulless barley in response to powdery mildew.
Plant defenses to biotrophic pathogens include an inherent defense system and an induced system that is activated after plants detect the attack. The former mainly includes cutin, waxy and lignin of plant cell walls, and small molecule disease-resistant compounds (e.g. fatty acids, phenolic substances, terpenoids and flavonoids). In the WGCNA of G7, the MEbrown module has the highest correlation, in which 350 genes were significantly enriched in phenylalanine metabolism, terpenoid backbone biosynthesis, sesquiterpenoid and triterpenoid biosynthesis, flavonoid biosynthesis, zeatin biosynthesis, isoflavonoid biosynthesis and so on, indicating that this set of genes might be responsible for the resistance of G7 in response to powdery mildew. For example, 12 genes were identified in the phenylalanine metabolic pathway; most of these genes were highly expressed at 0 hpi and 168 hpi in both G7 and Z13, and these genes had similar expression patterns in both cultivars (Fig. S4). These results showed that the plants might trigger similar gene regulation mechanisms and response processes to cope with the infection of powdery mildew. Therefore, it can be proposed that plants employ the same set of genes in defense pathways during the process of this plant-pathogen interaction and that the sequentially differential expression of these genes contributes to the defense performance of the plant-pathogen interactions.
The plant cell wall is impermeable to water and pathogens and is dramatically affected by cutin, suberine and wax deposition 32 . The cutin-wax barrier plays a crucial role in the interaction and adaptation of plants to various stresses, particularly by providing the first barrier to many pathogens 33 . In the MEblue module of G7, four genes were identified to be involved in cutin, suberine and wax biosynthesis (HVUL0H23092.2; HVUL4H38244.2; HVUL4H27323.2; HVUL1H15328.2, Fig. S5) encoding fatty acid omega-hydroxylase, omega-hydroxypalmitate O-feruloyl transferase, peroxygenase and aldehyde decarbonylase, respectively. In G7, these genes were expressed at a low level at 0 hpi, were upregulated at 6, 36, and 72 hpi, and returned to the 0 hpi level at 168 h. However, in Z13, three genes were expressed at very low levels and were not differentially expressed throughout the process of interaction with Bgh. Interestingly, further metabolomics analysis showed that the primary substrates of cutin, suberine and wax biosynthesis, lipids including glycerophospholipids and fatty acids accumulated more in G7 at all time points. These results might suggest that these genes and metabolites play important roles in the complete resistance of G7 to powdery mildew.
Phytohormones play crucial roles in a complex regulatory network that is essential for pathogen-induced responses. Interestingly, our results showed that the content of jasmonic acid (JA) was much higher in G7 than that in Z13 at 0 hpi, and the expression of genes associated with JA biosynthesis was also further examined. The data showed that several JA related genes, such as OPR, ACX, MPF2 and FADI, were more highly expressed in G7 at 0 hpi than they were in Z13 (Fig. 6A). Many studies have reported that JA plays an important role in defense responses against biotrophic pathogens 34 ; our findings are consistent with the conclusion that the accumulation of JA can have systematic effects on the resistance of powdery mildew in plants 22 .
Previous studies have suggested that the expression levels of PR are mediated through pathogen-induced signal-transduction pathways that are well regulated by phytohormones such as JA and MeJA 35 . The expression levels of 17 PRs were further analyzed in G7 and Z13 during the Qingke-Bgh interaction. Six genes (HVUL0H09210.2, HVUL3H06540.2, HVUL3H36018.2, HVUL5H11228.2, HVUL5H45926.2, HVUL7H00882.2) were induced by powdery mildew in G7 (Fig. S6), and the expression level of the PRs reached the maximum at 6 to 36 hpi. In the susceptive cultivar Z13, all 17 PRs were induced by powdery mildew (Fig. S6), and the expression level of the PRs reached the maximum at 168 hpi. The degree to which each PR gene was induced by Bgh differed in the two varieties. The result showed that the induced expression of PR genes was significantly different in G7 and Z13. The induced expression of PRs play more important roles in the resistant variety than in the susceptible one. These results indicated that powdery mildew infection had different effects on the 17 PRs, and the increased fold-changes in expression and the interval (response time) from inoculation to expression maxima were also different.
In conclusion, the combined transcriptome and metabolome analyses generated a set of data reflecting the dynamic defense responses of Tibetan hulless barley to powdery mildew. The defense responses involved primary metabolisms such as amino acid metabolism, carbohydrate metabolism, and lipids metabolism, and secondary metabolisms, including the biosynthesis of phenylpropanoids, flavone and flavonoids, phenolamides, and phytohormones. The gene expression and metabolic networks identified in this study provide new insights into the mechanisms of induced defense response in Qingke. The current findings will greatly improve our understanding of Qingke-Bgh interaction and provide clues for the development of resistant Tibetan hulless barley varieties.

Methods
Plant materials, Bgh inoculation and resistance. Two cultivated hulless barley G7 and Z13 were used in this study. G7 is a variety with complete Bgh resistance (showed complete penetration resistance to powdery mildew), while Z13 is sensitive to Bgh. In G7, after inoculated with Bgh, no significant changes in phenotype were observed during the long-time (168 h) process of Qingke-powdery mildew interaction. In Z13, the symptom of the disease is as follows: at 6 hpi, there was no significant changes; then at 36 hpi, there are tiny spots of fungal growth on the upper surface of the leaves. At 72 hpi, powdery mildew can be observed as little grayish spots and the spots become circular or irregular spots, which diameter reaches 1.0 mm. At 168 hpi, the patches connected into blocks and the spore density on the leaves were very thick with severe necrotic spots.
The hulless barley seedling cultivation and all experiments were carried out at the Tibet Academy of Agricultural and Animal Husbandry Sciences, Lhasa (Tibet, China). In the experiment, Qingke seedlings were grown in plant growth chambers under LD (18 h light:6 h dark) cycles (at 16-18 °C and a relative humidity of 80%). At the two-leaf and one needle stage, 30 well-developed seedlings were inoculated with Bgh. The growth temperature was adjusted to 20 °C after inoculation. The reactions to Bgh inoculation were observed visually based on the presence (or absence) of powdery mildew colonies and necrotic spots on the leaf surface. Subsequently, leaves from each individual plant were collected at 0, 6, 36, 72 and 168 hpi. sequenced on Illumina HiSeq X Ten platform. High-quality reads were generated for each library. Genome mapping of high-quality reads was carried out using the software TopHat2. The Tibetan hulless barley draft genome assembly (http://show.genebang.com/project/download?n=barley) was used as a reference for genome mapping of the 30 Qingke libraries.

Calculation of gene expression level and identification of differentially expressed genes.
The expression levels of Tibetan hulless barley genes were measured using the TPM (transcripts per million) method, and low-expression genes with TPM <1 were filtered in all 30 samples by the Kallisto program 36 . The reference used was the coding sequence of 36,151 Tibetan hulless barley genes. Differentially expressed gene (DEG) analysis between samples was performed using NOISeq 37 . The probability of 0.8 and [log 2 (Fold change)] greater than 1 were set as the thresholds for significant differential expression.
Gene expression pattern analysis was performed by Short Time-series Expression Miner software (STEM) 38 on the OmicShare tools platform (www.omicshare.com/tools). The parameters were set as follows: (1) Maximum Unit Change in model profiles between time points is 1; (2) Maximum output profile number is 20 (similar profiles will be merged); and (3) Minimum ratio of fold change of DEGs is no less than two.
WGCNA network analysis. The WGCNA R package was used to build GCNs, examine main network properties (hubs and modules), calculate the significance values of genes, and investigate the correlations between modules and powdery mildew resistance responses. The transcriptome datasets were filtered to remove any genes whose TMP value < 1 at any time point to remove these genes that introduce noise into the network analysis. Log 2 normalized TMP values were used to construct the coexpression networks using the WGCNA package 17 in R. Independent signed networks were generated from each Tibetan hulless barley cultivar G7 or Z13. The process of WGCNA network analysis was according to https://labs.genetics.ucla.edu/horvath/CoexpressionNetwork/ Rpackages/WGCNA/.
Widely Targeted Metabolome detection and data analysis. The freeze-dried leaf was ground using a mixer mill (MM 400, Retsch) with zirconia beads for 1.5 min at 30 Hz. 100 mg powder was extracted overnight at 4 °C with 1.0 ml 70% aqueous methanol. After centrifugation at 10, 000 g for 10 min, the extracts were absorbed and filtrated before LC-MS analysis.
The sample extracts were analyzed using an LC-ESI-MS/MS system (HPLC, Shim-pack UFLC SHIMADZU CBM30A system; MS, Applied Biosystems 4500 Q TRAP). The quantification of metabolites was accomplished using multiple reaction monitoring (MRM) analysis of a triple quadrupole-linear ion trap mass spectrometer (Q TRAP), API 4500 Q TRAP LC/MS/MS System. The analytical conditions and detailed operation parameters were as published before 39 .
For statistical analysis, the relative abundances of each metabolite were log transformed before analysis to meet normality. Dunnett's test was used to compare the abundance of each metabolite between different time points. Statistical analyses were performed using the SPSS 22.0 software package (IBM SPSS, Somers, NY, USA). Differentially expressed metabolites were identified using a combination of fold change and the VIP values of the OPLS-DA model. The VIP values ≥1 and [log 2 (Fold change)] greater than 1 were set as the threshold.