Bozeng Tang1, Caiyun Liu2, Zhiqiang Li1, Xixi Zhang1, Shaoqun Zhou3, Guo-Liang Wang4, Xiao-Lin Chen2, Wende Liu1. 1. State Key Laboratory for Biology of Plant Diseases and Insect Pests, Institute of Plant Protection, Chinese Academy of Agricultural Sciences, Beijing, China. 2. State Key Laboratory of Agricultural Microbiology and Provincial Hubei Key Laboratory of Plant Pathology, College of Plant Science and Technology, Huazhong Agricultural University, Wuhan, China. 3. Shenzhen Branch, Guangdong Laboratory of Lingnan Modern Agriculture, Genome Analysis Laboratory of the Ministry of Agriculture and Rural Affairs, Agricultural Genomics Institute at Shenzhen, Chinese Academy of Agricultural Sciences, Shenzhen, China. 4. Department of Plant Pathology, The Ohio State University, Columbus, OH, USA.
Abstract
Upon fungal and bacterial pathogen attack, plants launch pattern-triggered immunity (PTI) by recognizing pathogen-associated molecular patterns (PAMPs) to defend against pathogens. Although PTI-mediated response has been widely studied, a systematic understanding of the reprogrammed cellular processes during PTI by multi-omics analysis is lacking. In this study, we generated metabolome, transcriptome, proteome, ubiquitome and acetylome data to investigate rice (Oryza sativa) PTI responses to two PAMPs, the fungi-derived chitin and the bacteria-derived flg22. Integrative multi-omics analysis uncovered convergence and divergence of rice responses to these PAMPs at multiple regulatory layers. Rice responded to chitin and flg22 in a similar manner at the transcriptome and proteome levels, but distinct at the metabolome level. We found that this was probably due to post-translational regulation including ubiquitination and acetylation, which reshaped gene expression by modulating enzymatic activities, and possibly led to distinct metabolite profiles. We constructed regulatory atlas of metabolic pathways, including the defence-related phenylpropanoid and flavonoid biosynthesis and linoleic acid derivative metabolism. The multi-level regulatory network generated in this study sets the foundation for in-depth mechanistic dissection of PTI in rice and potentially in other related poaceous crop species.
Upon fungal and bacterial pathogen attack, plants launch pattern-triggered immunity (PTI) by recognizing pathogen-associated molecular patterns (PAMPs) to defend against pathogens. Although PTI-mediated response has been widely studied, a systematic understanding of the reprogrammed cellular processes during PTI by multi-omics analysis is lacking. In this study, we generated metabolome, transcriptome, proteome, ubiquitome and acetylome data to investigate rice (Oryza sativa) PTI responses to two PAMPs, the fungi-derived chitin and the bacteria-derived flg22. Integrative multi-omics analysis uncovered convergence and divergence of rice responses to these PAMPs at multiple regulatory layers. Rice responded to chitin and flg22 in a similar manner at the transcriptome and proteome levels, but distinct at the metabolome level. We found that this was probably due to post-translational regulation including ubiquitination and acetylation, which reshaped gene expression by modulating enzymatic activities, and possibly led to distinct metabolite profiles. We constructed regulatory atlas of metabolic pathways, including the defence-related phenylpropanoid and flavonoid biosynthesis and linoleic acid derivative metabolism. The multi-level regulatory network generated in this study sets the foundation for in-depth mechanistic dissection of PTI in rice and potentially in other related poaceous crop species.
Plants have developed elaborate defence mechanisms to prevent pathogen infection. Plants use plasma membrane‐located pattern recognition receptors (PRRs) to recognize pathogen‐derived pathogen‐associated molecular patterns (PAMPs; Jones and Dangl, 2006). This recognition is transduced to downstream signalling components to trigger a suite of molecular and physiological responses, collectively known as the pattern‐triggered immunity (PTI; Jones and Dangl, 2006). These responses include activation of mitogen‐activated protein kinase (MAPK) cascades in the cytosol and up‐regulation of transcription of defence‐related genes in the nucleus. To overcome the PTI‐mediated defence, some pathogens could secrete a number of effector proteins into the host cells. In turn, plants have evolved another layer of defence referred to as the effector‐triggered immunity (ETI), in which the plant cell cytoplasmic immune receptors recognize pathogen‐derived effectors (Cui et al., 2015; Tsuda and Katagiri, 2010).In plants, PTI serves as the first layer of defence conferring resistance to pathogens (Bigeard et al., 2015). PTI‐mediated defence is induced by PAMPs, including fungal‐derived chitin and bacterial‐derived flg22 (Boller and Felix, 2009). These PAMPs are recognized by specific receptors, such as the chitin receptor CERK1 and the flg22 receptor FLS2. PAMPs receptors recognition systems have been extensively studied in Arabidopsis and rice (Oryza sativa; Zipfel, 2008; Zipfel, 2014; Liu et al., 2014), which makes these species valuable scientific models to investigate the molecular mechanisms and cellular processes of PTI.Although PAMP receptors are specific, subsequent immune responses such as gene expression patterns are conserved (Bjornson et al., 2021; Boller and Felix, 2009; Jones and Dangl, 2006). In Arabidopsis, PTI‐mediated defence involves production of reactive oxygen species (ROS), transcriptional activation of pathogenesis‐related (PR) gene expression and secondary metabolites accumulation such as callose and phytoalexins (Boller and Felix, 2009). Excessive ROS is harmful to plants, and hence, plants produce various ROS detoxifying enzymes (Frederickson Matika and Loake, 2014). PAMP treatment induces PR genes including chitinase genes and biosynthetic genes of defensive secondary metabolites, such as flavonoids and lignin precursors (van Loon et al., 2006). During PTI‐mediated defence response, the plant hormones‐mediated signalling pathways, especially salicylic acid (SA), jasmonic acid (JA) and ethylene (ET) signalling, are also activated to mediate various cellular responses (Aerts et al., 2021). Like in Arabidopsis, similar PTI‐mediated defence occurs in rice against fungal and bacterial infection, such as ROS production, callose deposition, MAPK signalling cascade activation and the hormone signalling pathways activation (Liu et al., 2013, 2014; Yang et al., 2013). Up to now, some single or combined omics studies, including transcriptomic, proteomics and metabolomics, have been reported to investigate the effects of plant responses to invading pathogens or challenged by different PAMPs (Bassal et al., 2020; Bjornson et al., 2021; Chen et al., 2018; Jeon et al., 2020; Lovelace et al., 2018; Nobori et al., 2018; Winkelmüller et al., 2021). However, an integrative study of multi‐omics data is still lacking. Post‐translational modification, such as ubiquitination and acetylation, is known to modulate protein abundance and function (Withers and Dong, 2017). Therefore, it is important to explore how PAMP‐induced defence modulates rice immunity via post‐translational modifications (PTMs) as additional regulatory mechanisms. Moreover, most of the systemic PTI omics analyses were generated from the dicot plant Arabidopsis, and the conclusions from these existing studies may not be extended to rice and other monocot species.In this study, we generated and integrated metabolome, transcriptome, proteome, ubiquitome and acetylome data to elucidate rice responses to fungal and bacterial PAMPs, chitin and flg22. We revealed multiple regulatory landscapes in rice during PTI. We identified a set of defence responses that were common and distinct between flg22 and chitin treatment. Furthermore, we also mapped regulatory events in linoleic acid and alpha‐linolenic acid metabolism, and phenylpropanoid and flavonoid biosynthesis. This study provides evidences of molecular events at multiple layers and could facilitate further mechanistic dissection of PTI processes in rice and possibly other monocot species.
Results
Multi‐omics strategy for defence responses induced by chitin and Flg22
Plants trigger temporally coordinated immune responses such as accumulation of ROS, phosphorylation of MAPKs and expression of PR genes during PTI (Boller and Felix, 2009). We investigated these processes in 1‐week‐old rice seedlings treated with chitin or flg22 as in a previous study (Chen et al., 2018). Leaf fragments from rice seedlings were separately immersed in chitin or flg22 solution, and the levels of ROS were measured using the luminol chemiluminescence assay at different time points (from 0 to 20 min every 10 s). We found that ROS burst occurred rapidly after treatment of chitin or flg22 with the peak at 5 min (Figure 1a). MAPK3/6 phosphorylation level was significantly increased with an evident elevation at 0.5 h (Figure 1b). We also detected that the transcription of the PR and defence‐related genes (KS4, NAC4, PAL1 and PAL2) were elevated at 3 h after chitin and flg22 treatment (Figure 1c). Taken together, these results indicate that our experimental system is appropriate to monitor physiological responses during rice PTI.
Figure 1
Framework of integrative multi‐omics network analysis for chitin‐ and flg22‐triggered immune responses in rice. (a) Measurement of oxidative burst in rice seedlings treated with water control, chitin or flg22. (b) MAPK activation in rice seedlings treated with water control, chitin or flg22. Anti‐p42/44 antibody was used for Western blotting. Loading of total proteins was evaluated by Coomassie brilliant blue (CBB) staining. (c) Quantitative real‐time PCR analysis of pathogenesis‐related genes in response to chitin or flg22. For each sample, the ubiquitin gene OsUBQ1 was used as internal control. Error bars present the ± SD of three biological replicates. (d) Multiple‐omics approaches design and datasets summary.
Framework of integrative multi‐omics network analysis for chitin‐ and flg22‐triggered immune responses in rice. (a) Measurement of oxidative burst in rice seedlings treated with water control, chitin or flg22. (b) MAPK activation in rice seedlings treated with water control, chitin or flg22. Anti‐p42/44 antibody was used for Western blotting. Loading of total proteins was evaluated by Coomassie brilliant blue (CBB) staining. (c) Quantitative real‐time PCR analysis of pathogenesis‐related genes in response to chitin or flg22. For each sample, the ubiquitin gene OsUBQ1 was used as internal control. Error bars present the ± SD of three biological replicates. (d) Multiple‐omics approaches design and datasets summary.To comprehensively understand and compare the effects of the two PAMPs, we performed experiments to generate multi‐omics data. One‐week‐old rice seedlings were treated with chitin or flg22 for 3 h and subjected to different omics approaches including transcriptome, proteome, acetylome and metabolome analysis. Briefly, we detected 41 098 transcripts derived from 37 780 protein encoding genes, 7377 proteins, 839 acetylated proteins and 13 903 mass features (Figure 1d). In addition, we identified 2882 ubiquitinated proteins from our previous study (Figure 1d; Chen et al., 2018).To determine changes induced by PAMP treatment in rice, we compared chitin and flg22 responses in these five datasets using log2‐fold changes compared to the water‐treated control samples (Figure 2a,b). We found that all five P‐values of paired t‐tests derived from five datasets were lower than 0.0001, suggesting significant effects on rice were induced by the two PAMPs. To further determine the similarity between rice responses to chitin and flg22, Pearson’s correlation analysis was performed. Rice metabolites showed the least correlation between chitin and flg22 (R
2 = 0.34). By comparison, correlation coefficients of proteome and acetylome were 0.49 and 0.42 respectively. Ubiquitome and transcriptome showed higher R values, 0.62 and 0.63 (Figure 2c), respectively, indicating that rice responded to chitin and flg22 in a relatively similar manner at transcriptional and ubiquitination level.
Figure 2
Rice response to chitin and flg22 severely in multiple‐tiered regulation. (a) Overview of altered metabolites, transcripts, proteins and post‐translational modifications (PTMs) as ubiquitinated and acetylated proteins. The violin plots show the fold changes values in chitin‐ and flg22‐treated samples comparing to control samples. (b) Dot plots distribution of fold changes values of chitin‐treated and flg22‐treated samples compared to the control samples. The lines indicate the patterns of variation between chitin‐treated and flg22‐treated samples for each variable. (c) Pearson correlation analysis to determine the similarity between chitin‐treated and flg22‐treated samples. The x‐axis indicates fold changes values obtained for chitin‐treated samples, and y‐axis indicates fold changes values obtained for flg22‐treated samples. The R values indicate the degree of correlation.
Rice response to chitin and flg22 severely in multiple‐tiered regulation. (a) Overview of altered metabolites, transcripts, proteins and post‐translational modifications (PTMs) as ubiquitinated and acetylated proteins. The violin plots show the fold changes values in chitin‐ and flg22‐treated samples comparing to control samples. (b) Dot plots distribution of fold changes values of chitin‐treated and flg22‐treated samples compared to the control samples. The lines indicate the patterns of variation between chitin‐treated and flg22‐treated samples for each variable. (c) Pearson correlation analysis to determine the similarity between chitin‐treated and flg22‐treated samples. The x‐axis indicates fold changes values obtained for chitin‐treated samples, and y‐axis indicates fold changes values obtained for flg22‐treated samples. The R values indicate the degree of correlation.
Metabolite signatures during PTI
Metabolic response during PTI has been extensively studied in Arabidopsis (Misra et al., 2016; Schenke et al., 2011), but it is not clearly known in rice. To determine how PAMPs treatments impact rice metabolism, we generated metabolome dataset by performing a non‐targeted liquid chromatography‐mass spectrometry (LC‐MS) analysis with 10 biological replicates per treatment. (Dataset S1; Figure S1a, b). Compared to the control group, we identified 606 up‐regulated and 182 down‐regulated mass features (FDR < 0.05; fold‐change >2) in chitin‐treated samples, and 325 up‐regulated and 346 down‐regulated mass features in flg22‐treated samples (Figure 3a; Dataset S1), suggesting a stronger effect of chitin than flg22 on rice metabolism. We then used MetaboAnalyst version 4.0 (Chong et al., 2018) to annotate putative metabolites and associated pathways induced by chitin and/or flg22 in rice. In total, 113 significantly altered compounds were characterized, 97 of which showed distinct regulation patterns between chitin and flg22 treatments (Dataset S2), suggesting distinct rice responses to chitin and flg22 at the metabolome level. We then generated a predictive network to uncover metabolite–metabolite interactions extracted from STITCH 5 (Szklarczyk et al., 2016), according to the known reactions in PubMed. This enabled us to illustrate the potential functional relationships between detected metabolites and their associated metabolites in this study (Figure 3b). In chitin‐treated samples, compounds associated with energy metabolism, nucleotide metabolism, amino acids, pyridoxal phosphate (PLP) and serotonin were influenced. Meanwhile, in flg22‐treated samples, compounds associated with hexacosanoic acid, pyruvic acid and citric acid were influenced (Figure 3b). To determine metabolic pathways associated with the altered metabolites, we performed enrichment analysis using quantitative analysis of annotated compounds (Dataset S3). In total, 15 and 9 pathways were statistically enriched (P‐adj < 0.05) in the chitin‐ and flg22‐treated samples, respectively, supporting the possibility that rice responded to chitin more drastically than to flg22 at the metabolome level (Figure 3c). Starch and sucrose metabolism, galactose metabolism and amino sugar and nucleotide sugar metabolism were enriched in both chitin and flg22 samples (Figure 3c). Pathways such as phenylpropanoid biosynthesis, ascorbate and aldarate metabolism and tryptophan biosynthesis were enriched only in chitin‐treated samples. In contrast, pathways of diterpenoid biosynthesis, biotin metabolism, cutin, suberine and wax biosynthesis; biosynthesis of unsaturated fatty acids and alpha‐linolenic acid metabolism were enriched only in flg22‐treated samples (Figure 3c). We also found compounds significantly suppressed by both PAMPs, that is, abscisic aldehyde and hydroxychlorophyllide (Figure 3d). In contrast, catechin and tetrahydrofuran were induced by both treatments (Figure 3e). These compounds commonly impacted by both PAMPs could be defined as core metabolites that are responsive to PAMP treatment.
Figure 3
In‐depth metabolomics analysis defined distinct patterns in rice treated with chitin or flg22. (a) Abundances of metabolites features were quantified by untargeted metabolomics from 10 biologically independent replicates for each group from rice seedling samples. The cloud plots demonstrated metabolite features deemed significant by two‐sided ANOVA in contrasts (P ≤ 0.05) to an additional fold change cut‐off (fold change > 2) detected from positive polarity. The x‐axis indicated retention time of features, and the y‐axis indicated m/z values. The size of the circles indicated the degree of fold changes. The colours indicate up‐regulated (green) or down‐regulated (red) features comparing to the control (ck). (b) Metabolite–metabolite interaction networks highlighted the potential functional relationships between the sets of annotated metabolites, and how the altered compounds in chitin‐treated (red) and flg22‐treated (blue) samples impact to neighbour metabolites in the subnetworks. The chemical–chemical association for the metabolites network was extracted from STITCH. (c) Heatmap shows hierarchical clustering of the metabolic pathways from metabolomics data. Heatmap demonstrated z‐scores obtained from enrichment analysis by using fold changes values and P‐values from all the features detected from metabolomics. The colours indicated the significance as ‐log10 (FDR) of enrichment for each metabolic pathway in rice. (d) Boxplot shows intensities of detected metabolites which are down‐regulated in chitin or flg22 treatment. Each dot represents intensity of corresponding compounds detected in the metabolomics. All the compounds shown here are significantly down‐regulated (P‐adj < 0.05). (e) Boxplot shows intensities of detected metabolites which are up‐regulated in chitin or flg22 treatment. Each dot represents intensity of corresponding compounds detected in the metabolomics. All the compounds shown here are significantly up‐regulated (P‐adj < 0.05).
In‐depth metabolomics analysis defined distinct patterns in rice treated with chitin or flg22. (a) Abundances of metabolites features were quantified by untargeted metabolomics from 10 biologically independent replicates for each group from rice seedling samples. The cloud plots demonstrated metabolite features deemed significant by two‐sided ANOVA in contrasts (P ≤ 0.05) to an additional fold change cut‐off (fold change > 2) detected from positive polarity. The x‐axis indicated retention time of features, and the y‐axis indicated m/z values. The size of the circles indicated the degree of fold changes. The colours indicate up‐regulated (green) or down‐regulated (red) features comparing to the control (ck). (b) Metabolite–metabolite interaction networks highlighted the potential functional relationships between the sets of annotated metabolites, and how the altered compounds in chitin‐treated (red) and flg22‐treated (blue) samples impact to neighbour metabolites in the subnetworks. The chemical–chemical association for the metabolites network was extracted from STITCH. (c) Heatmap shows hierarchical clustering of the metabolic pathways from metabolomics data. Heatmap demonstrated z‐scores obtained from enrichment analysis by using fold changes values and P‐values from all the features detected from metabolomics. The colours indicated the significance as ‐log10 (FDR) of enrichment for each metabolic pathway in rice. (d) Boxplot shows intensities of detected metabolites which are down‐regulated in chitin or flg22 treatment. Each dot represents intensity of corresponding compounds detected in the metabolomics. All the compounds shown here are significantly down‐regulated (P‐adj < 0.05). (e) Boxplot shows intensities of detected metabolites which are up‐regulated in chitin or flg22 treatment. Each dot represents intensity of corresponding compounds detected in the metabolomics. All the compounds shown here are significantly up‐regulated (P‐adj < 0.05).
Transcriptomic signature during PTI
We reasoned that the metabolic changes in PAMPs‐treated rice seedlings were underpinned by alteration in gene expression. Therefore, we compared transcriptomes of rice seedlings treated by chitin or flg22. After alignment of the clean reads to the rice reference genome (IRGSP‐1.0), the read counts were used to validate data quality by Pearson correlation and principal component analysis (PCA; Figure S2a,b). For differentially expressed genes analysis (DEGs) (FC > ±2, P‐adj < 0.05), we identified 531 genes that were specifically up‐regulated by chitin, 656 genes that were up‐regulated only by flg22 and 287 genes that were up‐regulated by both PAMPs. We also identified 240 down‐regulated genes in response to chitin specifically, 260 down‐regulated genes in response to flg22 specifically and 120 down‐regulated genes in response to both PAMPs (Figure 4a,b; Dataset S4). Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis (P‐adj < 0.05) showed that pathways associated with glycolysis/gluconeogenesis, phenylpropanoid biosynthesis, flavonoid biosynthesis, fructose and mannose metabolism and carbon metabolism were induced by both PAMPs (Figure 4c). Alanine, aspartate and glutamate metabolism and pentose phosphate pathway were specifically responsive to chitin (Figure 4c,d). Linoleic acid metabolism and alpha‐linoleic acid metabolism were only significantly induced by flg22 but not chitin, which was consistent with the results of metabolome analysis (Figure 3c). Gene ontology terms significantly enriched among DEGs (FDR < 0.05) were also identified (Figure S2c,d).
Figure 4
The rice seedlings transcriptome is substantially altered upon chitin or flg22 treatment. (a) Volcano plot comprises of chitin treatment to ck, and flg22 treatment to ck from RNA‐seq analysis. The colours indicated genes whether they significantly differentially expressed genes under the cut‐off of (P < 0.05) or/and fold changes (FC > 2). (b) Venn diagram shows the numbers of differentially expressed genes identified in chitin‐treated and flg22‐treated samples. The numbers indicated the total numbers of either up‐regulated or down‐regulated genes. The overlapping differentially expressed genes were also indicated. DEGs are defined under the cut‐off of P‐adj < 0.05 and FC > 2. (c) Network plots the detailed survey to visualize detailed enrichment of metabolic pathways in rice impact by DEGs in chitin‐treated and flg22‐treated samples. The colours indicate the comparison groups derived from RNA‐seq analysis. (d) Dot plot shows the P‐values as significances of enrichment of the pathways enriched by DEGs. The size of dots indicates the gene ratio of DEGs in each pathway. The colour indicates the degrees of significance as P‐adj values. Only pathways showing P‐adj less than 0.05 are visualized and analysed. (e) Dot plot shows fold changes of DEGs in RNA‐seq of chitin‐treated and flg22‐treated samples compared to control samples. The lines indicated the paired values detected in chitin‐treated and flg22‐treated samples.
The rice seedlings transcriptome is substantially altered upon chitin or flg22 treatment. (a) Volcano plot comprises of chitin treatment to ck, and flg22 treatment to ck from RNA‐seq analysis. The colours indicated genes whether they significantly differentially expressed genes under the cut‐off of (P < 0.05) or/and fold changes (FC > 2). (b) Venn diagram shows the numbers of differentially expressed genes identified in chitin‐treated and flg22‐treated samples. The numbers indicated the total numbers of either up‐regulated or down‐regulated genes. The overlapping differentially expressed genes were also indicated. DEGs are defined under the cut‐off of P‐adj < 0.05 and FC > 2. (c) Network plots the detailed survey to visualize detailed enrichment of metabolic pathways in rice impact by DEGs in chitin‐treated and flg22‐treated samples. The colours indicate the comparison groups derived from RNA‐seq analysis. (d) Dot plot shows the P‐values as significances of enrichment of the pathways enriched by DEGs. The size of dots indicates the gene ratio of DEGs in each pathway. The colour indicates the degrees of significance as P‐adj values. Only pathways showing P‐adj less than 0.05 are visualized and analysed. (e) Dot plot shows fold changes of DEGs in RNA‐seq of chitin‐treated and flg22‐treated samples compared to control samples. The lines indicated the paired values detected in chitin‐treated and flg22‐treated samples.During PTI, the host recognizes PAMPs by receptors located on the surface of plant cells, and transmits the immune signal through receptor‐like kinase (RLK), MAPK cascade and CDPK (calcium‐dependent protein kinase), activating the expression of transcription factors and subsequent defence‐related genes, such as hormone signalling pathway genes and secondary metabolism pathway genes (Boller and Felix, 2009; Jones and Dangl, 2006). Recently, a set of more than 1,000 genes in Arabidopsis thaliana has been reported to be the core PAMP‐responsive genes at transcription level (Bjornson et al., 2021). For a comparison analysis, we extracted the orthologous genes between A. thaliana and rice, and compared our DEGs to the orthologues of the Arabidopsis core PAMP‐responsive genes. A total of 51 genes shared between rice and Arabidopsis were identified as they were differentially expressed in either chitin‐ or flg22‐treated samples in our study (Dataset S4). Several of these genes encode WRKY domain‐containing proteins, phosphatases for MAPK and Ca2+ modulation proteins, indicating they are key factors in PTI in both species. To shed light on the underlying response to PAMPs at transcriptional level in rice, we extracted the DEGs associated with defence response, and found that most of which were similar in responding to chitin and flg22 (Figure 4e; Figure S3a; Dataset S5). Many RLK genes, including receptor‐like cytoplasmic kinases (RLCKs) and wall‐associated kinase (WAK) genes, were induced by both PAMPs. Calcium signalling pathway genes, such as some calmodulin‐like protein (CML) genes, were also induced by both PAMPs. Additional to that, 24 WRKY TFs were detected in the transcriptome, nearly all of which were induced by both PAMPs, including WRKY62 and WRKY76, which were known to play key roles in rice blast disease resistance (Liu et al., 2016). Other PAMPs‐induced TF genes include homeobox TFs, bZIP TFs, NAC TFs, Dof family TFs, heat stress TFs, BSRD1, JAMYB and CEF1/OsMYB103.The salicylic acid (SA), jasmonic acid (JA) and ethylene (ET) signalling pathways are three key players in plant PTI response (Berens et al., 2017; De Vleesschauwer, et al., 2014). We found many ET signalling pathway genes (such as ERFs, DREB1s, PLTs and AP2/EREBP89) and JA signalling pathway genes (such as OsDR10, JAMYB, AOS2 and OsOPR7) that were up‐regulated by both PAMPs (Figure 4e; Figure S3a; Dataset S5). Other hormone pathways, including IAA, ABA, GA and BR signalling pathways, may also affect SA‐JA‐ET pathways to regulate plant defence response (Figure 4e; Figure S3a; Dataset S5) (De Vleesschauwer, et al., 2014). Similarly, many genes involved in these hormone pathways were induced by both PAMPs (Figure 4e; Figure S3b; Dataset S5). We also found many redox reaction‐related genes were induced by both PAMPs, whereas many genes encoding laccases, glutathiones, catalases, oxidase and peroxiredoxin were repressed by both PAMPs. However, among the identified PR genes in transcriptome, PR1‐73, CHT11/PR3 and CHT6/PR‐3 were induced by chitin, PR10B, PR1‐73, PR5, CHT11/PR3, CHT6/PR‐3 and CHT8 were induced by flg22 (Figure S3a), indicating difference between chitin‐ and flg22‐induced PR gene expression.Taken together, the transcriptome analysis suggests that rice response to chitin and flg22 is in a similar manner to regulate gene expression, because many defence‐related genes are commonly induced in response to both PAMPs, which serves as core factors tightly associated with defence‐related pathways. Furthermore, numerous genes in several metabolic pathways such as glycolysis/gluconeogenesis, phenylpropanoid biosynthesis and flavonoid biosynthesis are significantly affected at transcriptional level (Figure S3c), and this is consistent with result of our comparative metabolomics analyses.
Integrative analysis of proteome and transcriptome
We employed iTRAQ approach to perform proteomics analyses for chitin, flg22 and water‐treated rice seedling samples (Figure S4a–e). In total, 7,387 detected proteins were used for data analysis (Dataset S6). We compared differentially expressed proteins (DEPs) between chitin and flg22 responses (FC > ±1.5, P‐adj < 0.05; Figure 5a,b; Dataset S7). As shown in Figure 5b, we found that 278, 226 and 76 proteins were significantly induced by chitin, flg22 or both respectively. In contrast, 297, 274 and 107 proteins were significantly suppressed by chitin, flg22 or both respectively. Further gene ontology (GO) enrichment analysis suggested that DEPs by both chitin and flg22 were associated with defence‐related processes, such as oxidation reduction, response to stress and response to stimulus (Figure 5c). Proteins associated with oxidoreductase activity, transition metal ion binding and cytoplasmic vesicle were shown to be enriched in responding to chitin specifically. Further KEGG metabolic pathway enrichment analysis suggested that phenylpropanoid biosynthesis, ether lipid metabolism, biosynthesis of amino acids, carbon fixation in photosynthetic organisms and inositol phosphate metabolism were enriched in both PAMP treatments (Figure 5d,e). Metabolic pathways of mRNA surveillance and arginine biosynthesis were only significantly enriched in chitin‐induced samples (P‐adj < 0.1) (Figure 5d,e). Metabolic pathway of fructose and mannose metabolism was significantly enriched only in flg22‐induced samples (P‐adj < 0.1) (Figure 5d,e). These results clearly suggested that, similar to the rice transcriptome, chitin and flg22 induced similar defence response at translational level.
Figure 5
Global protein profiling revealed a strong influence of PAMPS treatment on the rice proteome. (a) Volcano plot shows both ‐log10(P‐adj) and log2‐fold changes of all the proteins detected from proteomics in chitin‐treated and flg22‐treated samples, after comparing to ck samples. The back dots indicated differentially expressed proteins (DEPs) significantly (P‐adj < 0.05, FC > 1.5). (b) Venn diagram shows numbers of differentially expressed proteins detected in chitin‐treated and flg22‐treated samples. (c) Heatmap shows the normalized ‐log10 (P‐adj) values to demonstrate enrichment analysis using DEPs and defined gene ontology terms enriched by DEPs. (d) Dot plots show significances of enrichment for the metabolic pathways detected by using DEPs. (e) Network plot showed detailed information of pathways enriched by DEPs.
Global protein profiling revealed a strong influence of PAMPS treatment on the rice proteome. (a) Volcano plot shows both ‐log10(P‐adj) and log2‐fold changes of all the proteins detected from proteomics in chitin‐treated and flg22‐treated samples, after comparing to ck samples. The back dots indicated differentially expressed proteins (DEPs) significantly (P‐adj < 0.05, FC > 1.5). (b) Venn diagram shows numbers of differentially expressed proteins detected in chitin‐treated and flg22‐treated samples. (c) Heatmap shows the normalized ‐log10 (P‐adj) values to demonstrate enrichment analysis using DEPs and defined gene ontology terms enriched by DEPs. (d) Dot plots show significances of enrichment for the metabolic pathways detected by using DEPs. (e) Network plot showed detailed information of pathways enriched by DEPs.Rice protein abundance changed by PAMP treatment could be reflected by changes at the transcript level. To test this idea, we performed correlation analysis by comparing the relative levels as log2FC values of transcripts and proteins affected by two PAMPs but found very weak correlation between transcriptome and proteome responses in both chitin and flg22 treatments (R2 = 0.15/0.21; Figure 6a). We detected a number of genes showing convergent pattern of mRNA and proteins levels as core responsive factors towards both PAMPs, and further enrichment analysis suggests that they are tightly associated with membrane‐bounded vesicle, catalytic activity and lipid metabolic process (Figure 6b). Taken together, these results suggest that there are significant post‐transcriptional and/or post‐translational events that modulate protein abundances during rice responses to chitin and flg22.
Figure 6
Integrative analysis of transcriptomics and proteomics. (a) mRNA/protein expression fold changes in chitin‐treated and flg22‐treated samples compared to control. mRNAs/proteins differentially expressed (|log2[FC]| > 1; P‐adj < 0.05). (b) mRNAs/proteins differentially expressed in consistent or inconsistent patterns are coloured as indicated. The enriched gene ontology function enriched by groups of mRNAs/proteins is shown for each category. List of GO terms enriched in the group of proteins that are significantly induced at both mRNA and protein levels, or proteins that are significantly suppressed only at the protein level. (c) Heatmap shows integrated enrichment analysis of metabolic pathways using DEGs, DEPs and result from metabolomics pathways studies, using ‐log10 (P‐values) from enrichment analysis derived from each analysis. Regulatory networks of metabolic pathways of rice impact by chitin (d) and flg22 (e) at protein level and transcript level. The up‐regulated or down‐regulated at both levels were indicated by colours.
Integrative analysis of transcriptomics and proteomics. (a) mRNA/protein expression fold changes in chitin‐treated and flg22‐treated samples compared to control. mRNAs/proteins differentially expressed (|log2[FC]| > 1; P‐adj < 0.05). (b) mRNAs/proteins differentially expressed in consistent or inconsistent patterns are coloured as indicated. The enriched gene ontology function enriched by groups of mRNAs/proteins is shown for each category. List of GO terms enriched in the group of proteins that are significantly induced at both mRNA and protein levels, or proteins that are significantly suppressed only at the protein level. (c) Heatmap shows integrated enrichment analysis of metabolic pathways using DEGs, DEPs and result from metabolomics pathways studies, using ‐log10 (P‐values) from enrichment analysis derived from each analysis. Regulatory networks of metabolic pathways of rice impact by chitin (d) and flg22 (e) at protein level and transcript level. The up‐regulated or down‐regulated at both levels were indicated by colours.
Convergence and divergence of rice metabolic pathway responses in chitin and Flg22 PTI at multiple regulatory layers
We observed distinct responses between chitin and flg22 PTI in rice metabolome, but similar responses in rice transcriptome and proteome. By combining transcriptome, proteome and metabolome results (Figure 6c), we identified a number of pathways enriched in both chitin and flg22 treatment, such as glycolysis/gluconeogenesis, phenylpropanoid biosynthesis and amino sugar and nucleotide sugar metabolism. Among them, phenylpropanoid biosynthesis pathway was induced in all samples, except in the metabolome of flg22‐treated samples. Flavonoid biosynthesis, which was downstream of the core phenylpropanoid biosynthesis, was induced in both PAMPs significantly at transcription level. Consistently, we also observed a number of metabolites altered in flavonoid biosynthesis (Dataset S3). Interestingly, flavonoid biosynthesis pathway was significantly enriched among DEGs, but not among DEPs or altered metabolites, suggesting that there may be additional post‐transcriptional regulation mechanism to modulate this process. Pathways associated with carbon metabolism and 2‐oxocarboxylic acid metabolism were enriched among DEGs and DEPs but not for metabolites, possible because less compounds were detected. Moreover, we found flg22‐specific responsive pathways, such as linoleic acid metabolism, alpha‐linolenic acid metabolism, biosynthesis of non‐saturated fatty acids, diterpenoid, sesquiterpenoid and triterpenoid biosynthesis, biotin metabolism and riboflavin metabolism, suggesting that they could play key roles in responses to flg22. We visualized enriched pathways by providing DEGs and DEPs identified (Figure 6d,e). We observed that the protein abundance and mRNA associated with phenylpropanoid biosynthesis were significantly changed when treated with chitin or flg22, implying that phenylpropanoid biosynthesis pathway was important in rice to respond to PAMPs. Strikingly, starch and sucrose metabolism were controlled by increase at both mRNA and protein levels upon chitin treatment (Figure 6d). These results indicated that chitin and flg22 induced similar defence response by modulating phenylpropanoid biosynthesis and starch and sucrose metabolism. The ether lipid metabolism was regulated by decreasing protein level upon chitin induction (Figure 6d). Inositol phosphate metabolism was influenced at both mRNA and protein levels upon flg22 induction (Figure 6e).Taken together, these results suggested that rice positively activated phenylpropanoid biosynthesis, flavonoid biosynthesis and glycolysis/gluconeogenesis in response to chitin and flg22 in a similar manner at transcriptional and translational level, but the responses at metabolic level were distinctive.
Dynamic post‐translational protein modulation during chitin and Flg22 PTI
Previously, we discovered dynamic protein ubiquitination in rice treated with chitin and flg22 in a comparative ubiquitome study (Chen et al., 2018), suggesting ubiquitination is important for protein abundance and function during PTI response. To test this idea, we processed our previous ubiquitome dataset (iProX, project ID IPX0001089000; Dataset S8) by imputing missing values and normalization of the counts of detected peptides (Figure S5a–e). We set up cut‐off (P < 0.05 and FC > ±1.5) to define differentially ubiquitinated proteins (DUPs) in chitin‐ and flg22‐treated samples (Dataset S9). We then integrated proteome and ubiquitome datasets, and observed that many DEPs detected in proteomics were also differentially ubiquitinated (Figure 7a,b). This result suggested that rice may regulate expression or translation of these proteins via ubiquitination, which may explain the inconsistent responses between transcriptome and proteome. Among them, 23 and 14 up‐regulated proteins were significantly more ubiquitinated in response to chitin and flg22 respectively. Meanwhile, 20 up‐regulated proteins were significantly less ubiquitinated in response to chitin compared to 9 proteins that demonstrated the same abundance and ubiquitination pattern in response to flg22 treatment. Furthermore, 10 down‐regulated proteins were significantly more ubiquitinated in chitin treatment compared to 16 such proteins in the flg22‐treated group (Figure 7a,b). This analysis captured a regulatory signature to protein expression in rice as the increased chitin‐responsive proteins could be potentially degraded via ubiquitination after 3 hpi, and decreased proteins maybe increased by lowered ubiquitination rather than up‐regulated transcription. Consistent with this assumption, these proteins were tightly involved in processes such as ribosome and translational elongation, which were associated with protein expression or translational regulation. To gain further insights into other PTMs during PTI in rice, we performed genome‐wide acetylome analysis. We deployed three biological replicates for chitin, flg22 and water‐treated rice seedling samples, following correlation analysis and normalization (Figure S6a‐e). In total, we found 10 572 acetylated sites corresponding to 839 protein‐encoding genes in rice (Dataset S10). Differentially acetylated proteins (DEAs) were identified using normalized intensities of acetylated proteins (Figure 7c‐e; Dataset S11). There were 34 and 110 proteins significantly more acetylated in chitin‐treated and flg22‐treated rice tissues respectively (P‐adj < 0.05, FC > ±2). There were 43 and 6 proteins significantly less acetylated when treated by chitin and flg22 respectively. Based on enrichment analysis, we observed that the differentially acetylated proteins in response to chitin and flg22 were predominantly involved in ribosomal subunit, ribonucleoprotein complex, protein‐containing complex, nucleosome, DNA packaging complex, translational elongation and translation (Figure S7a,b). These results clearly suggested that the acetylation process modulated proteins that played important roles in translational regulation to regulate gene expression. We also found a number of histone proteins in differentially acetylated proteins (Figure 7f), suggesting that epigenetics regulation could also be involved in PTI response in rice. These results highlighted sophisticated post‐translational modification processes are deployed to regulate protein abundance in rice in response to PAMPs.
Figure 7
Correlation between protein and ubiquitination or acetylation. Proteins/ubiquitinated proteins expression fold changes in chitin treated (a) and flg22 treated (b) were demonstrated after comparing to control samples. Proteins/ubiquitinated proteins differentially expressed (|log2[FC]| > 0.83; P‐adj < 0.05) in both, either and neither (unchanged) ubiquitome and proteome studies were grouped and colour coded. The numbers of proteins/proteins ubiquitinated are shown for each category. List of enriched GO terms in the box are significantly affected proteins or ubiquitinated proteins by chitin or flg22 treatment. Green dots represent genes showing up‐regulated in proteomics and more ubiquitinated. Red dots represent genes showing down‐regulated in proteomics and less ubiquitinated. Blue dots represent genes showing up‐regulated in proteomics but less ubiquitinated. Orange dots represent genes showing down‐regulated in proteomics but more ubiquitinated. Volcano plot shows the differentially acetylated proteins detected in chitin‐treated (c) and flg22‐treated (d) samples. (e) Venn diagram shows the numbers of differentially acetylated proteins detected in acetylome by using cut‐off as P‐adj < 0.05 and FC > 2. (f) The heatmaps show relative acetylation levels of histone proteins, transcription factors and ribosome proteins in rice responding to chitin and flg22 after comparing to control samples.
Correlation between protein and ubiquitination or acetylation. Proteins/ubiquitinated proteins expression fold changes in chitin treated (a) and flg22 treated (b) were demonstrated after comparing to control samples. Proteins/ubiquitinated proteins differentially expressed (|log2[FC]| > 0.83; P‐adj < 0.05) in both, either and neither (unchanged) ubiquitome and proteome studies were grouped and colour coded. The numbers of proteins/proteins ubiquitinated are shown for each category. List of enriched GO terms in the box are significantly affected proteins or ubiquitinated proteins by chitin or flg22 treatment. Green dots represent genes showing up‐regulated in proteomics and more ubiquitinated. Red dots represent genes showing down‐regulated in proteomics and less ubiquitinated. Blue dots represent genes showing up‐regulated in proteomics but less ubiquitinated. Orange dots represent genes showing down‐regulated in proteomics but more ubiquitinated. Volcano plot shows the differentially acetylated proteins detected in chitin‐treated (c) and flg22‐treated (d) samples. (e) Venn diagram shows the numbers of differentially acetylated proteins detected in acetylome by using cut‐off as P‐adj < 0.05 and FC > 2. (f) The heatmaps show relative acetylation levels of histone proteins, transcription factors and ribosome proteins in rice responding to chitin and flg22 after comparing to control samples.
Multi‐omics analysis uncovered detailed convergence and divergence in metabolic pathways associated with defence response
Our omics analysis enabled us to understand the underlying connections between the metabolites and its upstream enzymatic reactions controlled by the gene expression or PTMs of proteins. We further focused on phenylpropanoid biosynthesis and flavonoid biosynthesis because they were defence‐associated and responsive to both PAMPs in our dataset (Dixon et al., 2002; Treutter, 2005). We also focused on the pathway glycolysis/gluconeogenesis as it showed the highest degree of significance in transcriptome response, and all the relevant compounds were detected in metabolome. In addition, we also selected the linoleic acid metabolism for further in‐depth analysis, as it showed specific response to flg22. Using simplified models of metabolic pathways from KEGG and Pathview (Luo et al., 2017), we mapped the variables with quantitative changes occurred in at least one dataset, in our analysis (Figure 8).
Figure 8
High‐resolution mapping of regulatory network to metabolic pathways. Simplified metabolic flow schemes described changes in metabolites, and enriched enzymes associated with transcripts, proteins, ubiquitinated proteins and acetylated proteins in phenylpropanoid biosynthesis (a), flavonoid biosynthesis (b), glycolysis/gluconeogenesis (c), linoleic acid metabolism (d) and alpha‐linolenic acid metabolism (e). The colours indicated rates of up‐ or down‐regulation of groups of enzymes.
High‐resolution mapping of regulatory network to metabolic pathways. Simplified metabolic flow schemes described changes in metabolites, and enriched enzymes associated with transcripts, proteins, ubiquitinated proteins and acetylated proteins in phenylpropanoid biosynthesis (a), flavonoid biosynthesis (b), glycolysis/gluconeogenesis (c), linoleic acid metabolism (d) and alpha‐linolenic acid metabolism (e). The colours indicated rates of up‐ or down‐regulation of groups of enzymes.Increased abundance of enzymes [6.2.1.12], [1.14.1311] in phenylpropanoid metabolism appeared to be mediated by increased transcription, which in turn, partially increased the abundance of the corresponding metabolites (Figure 8a). We observed consistent regulation in response to both PAMPs at transcriptional and PTMs levels, but divergent pattern at the metabolic level, as eight compounds were differentially altered between chitin and flg22 PTI. Seven co‐induced metabolites such as phenylalanine, p‐cinnamic acid, 5‐hydroxyferulic acid and secoisolariciresinol were associated with significant increases in mRNA but reduced ubiquitination levels, implying a coordinated regulation. We found genes associated with enzyme [1.2.1.44] were translationally down‐regulated by both PAMPs, which correlated with increased abundance of caffeylaldehyde and coniferylaldehyde, but decreased cinnamaldehyde and sinapaldehyde. We did not detect significant changes in cinnamoyl‐CoA, and its derivatives, p‐cinnamoyl‐CoA, caffeoyl‐CoA and feruloyl‐CoA (Figure 8a). However, in flavonoid biosynthesis, genes associated with enzymatic reaction [2.3.1.74], a vital step to process downstream compounds, were down‐regulated by both PAMPs. And this led to significantly increased metabolites, such as pinocembrin chalcone, naringenin chalcone and eriodictyol (Figure 8b).In glycolysis/gluconeogenesis pathway, alpha‐D‐glucose‐1P, alpha‐D‐glucose‐6P, oxaloacetate and pyruvate were significantly up‐regulated in chitin‐treated tissues. Meanwhile, D‐glucose, glyceraldehyde‐3P, glycerone‐P and acetyl‐CoA were down‐regulated by both PAMPs. The abundance of alpha‐D‐glucose‐1P in rice was different between chitin and flg22‐treated samples. However, its derivative, beta‐D‐Fructose, was increased via enzymatic reaction [2.7.1.11], which is associated with elevated mRNA level. Moreover, reduced ubiquitination levels were detected in flg22‐responsive proteins (Figure 8c). Likewise, we observed lower abundance of flg22‐responsive oxaloacetate than that in chitin, and rice modulates genes associated with enzymatic reaction [41.1.49] to increase flg22‐responsive pyruvate to the same level as chitin (Figure 8c).Interestingly, linolenic acid‐associated pathways were specifically activated in flg22‐treated samples. During flg22‐responsive process, gene expression was negatively regulating abundance of metabolites (Figure 8d). By contrast, in alpha‐linolenic acid metabolism, almost all enzymatic reactions were positively activated by enhanced transcription in flg22‐treated samples and this correlated with increased abundance of eight compounds in this pathway (Figure 8e).Taken together, our integrated analysis provides dynamic metabolic landscape controlled by multilayers including transcription, translation and post‐translational protein modification that serves as a foundation for in‐depth molecular characterization by researchers.
Discussion
Previous studies have applied single or dual omics to uncover the molecular mechanisms involved in plant–pathogen interactions (Bassal et al., 2020; Bjornson et al., 2021; Chen et al., 2018; Jeon et al., 2020; Lovelace et al., 2018; Nobori et al., 2018, 2020; Winkelmüller et al., 2021). In this study, we performed a multi‐omics analysis consisting of five omics datasets including metabolome, transcriptome, proteome, ubiquitome and acetylome to systematically elucidate PTI triggered by two different PAMPs in rice. The analysis summarized a landscape of detailed convergence and divergence between chitin and flg22 responses in rice. Compared to previous single‐ or dual‐omics analyses, our integrated multi‐omics analysis revealed dynamic and distinct multi‐layered regulation during chitin and flg22 PTI and set the foundation for further characterization in depth. We found chitin‐triggered PTI responses were more similar to flg22 at transcriptional and translational levels, but much stronger and different from flg22 at PTMs and metabolome levels, which were overlooked by previous analyses. We also deduced detailed regulatory mechanisms that may cause differences in several plant defence‐related metabolic pathways in response to the two PAMPs. Moreover, the linoleic acid‐associated pathway was identified as an flg22‐specific responsive pathway through this multiple‐omics analysis.We have clearly shown that rice transcriptomic responses to flg22 and chitin in a very similar manner (Figures 2 and 4), which was in line with a recent study in Arabidopsis thaliana that found a suite of conserved response to flg22 and chitin (Bjornson et al., 2021). In this study, we summarized details of PAMPs‐induced genes in rice and correlated them with defence functions. A considerable number of genes, especially those involved in early defence signalling pathways, were commonly induced by both PAMPs at the transcription level (Figure 4). Furthermore, our RNA‐seq analysis clearly provided evidence of how different PAMPs affected rice metabolic pathways in a convergent manner, particularly for glycolysis/gluconeogenesis, phenylpropanoid biosynthesis and flavonoid biosynthesis. Consistent with the metabolomic analysis, we uncovered that linoleic acid‐associated pathways in rice were very likely to be specifically affected by flg22. Additional to that, many defence‐related genes were found among chitin‐ and flg22‐induced DEGs, including genes encoding RLKs, MAPKs and CDPKs, transcription factors, hormone signalling components and secondary metabolic genes (Figure 4e; Figure S3). Interestingly, some key growth–defence balance master regulator genes, including OsALDH2B1 (Ke et al., 2020), BSRD1 (Li et al., 2017), GRF9 (Zhang et al., 2018), HD16 (Nemoto et al., 2018), DWARF11 (Zhou et al., 2017), GLW7/SPL13 (Si et al., 2016) and EBR1 (E3 ubiquitin ligase interacts with OsBAG4 for balance; You et al., 2016), were identified as DEGs (Figure 4e; Figure S8a,b). These findings highlighted a genome‐wide landscape of growth–defence trade‐off balance during rice PTI response.We showed that rice metabolic responses to chitin and flg22 were distinct, as there were more pathways significantly enriched in chitin‐treated samples than that in flg22‐treated ones. This was consistent with our phenotype assay, from which we observed higher level of ROS in flg22‐treated samples than that in chitin‐treated ones, and higher degree of phosphorylation of MAPK in chitin‐treated samples than that in flg22‐treated ones (Figure 1). This was probably due to a more adaptive response occurring during plant responses to chitin. By further predictive analysis, we found pyridoxal phosphate (PLP) and serotonin were highly affected in chitin‐treated samples. PLP is an active form of vitamin B6, which was reported to modulate cellular antioxidant capacity for plant immune system (Zhang et al., 2014a, 2014b). Serotonin is a primary amino compound that is the 5‐hydroxy derivative of tryptamine, which involves in plant innate immunity by affecting cellular redox status and callose accumulation, the latter functions as a mechanical barrier against pathogens (Nehela and Killiny, 2020; Zhao et al., 2015). We also defined a group of novel metabolites discovered in the metabolome as putative PAMPs‐responsive compounds, as they showed convergent pattern to both PAMPs, such as abscisic aldehyde, octadcatrienoyl‐CoA and catechin (Figure 3). Abscisic acid (ABA) is a plant hormone involved in seed development and responds to various environmental stresses. Oxidation of abscisic aldehyde is the last step of ABA biosynthesis and is catalysed by aldehyde oxidase (Seo et al., 2000). These findings highlighted the signatures of metabolites in rice responding to chitin and flg22, and explained the more severe response to chitin at metabolomic level.Further integrative analysis of transcriptome with proteome in rice uncovered an inconsistent pattern between mRNA and protein abundance in both PAMPs treatments, although we concluded that rice response to both PAMPs was similar at the translational level. One of the explanations to this inconsistency was that RNA‐seq could detect mRNA at higher sensitivity than protein detection with the iTRAQ‐based proteomics platform. However, further analysis using the previous ubiquitome dataset (Chen et al., 2018) suggested that a protein turnover regulation may play a key role in this disparity, as we observed that a considerable number of DEPs were differentially ubiquitinated (Figure 7). The differentially ubiquitinated proteins were significantly enriched in ribosome functions and translation‐related regulations, indicating their potential roles in further defence response after 3 h treatment of chitin and flg22. We characterized 180 differentially acetylated proteins and found more DEAs in chitin‐treated samples than that in flg22‐treated samples. By contrast, further analysis with acetylome suggested a convergent pattern of core DEAs between chitin‐ and flg22‐treated samples, including histone‐related genes, transcription factor‐encoding genes and ribosome‐encoding genes (Figure 7). Further gene ontology enrichment analyses confirmed the importance of these DEAs as core units acting as gene expression regulators to modulate response to PAMPs in rice (Figure S7). Acetylome analysis not only summarized distinct responses between chitin‐ and flg22‐treated samples, in general, but also uncovered the core units that respond to both PAMPs. PTMs analysis demonstrated that there is a sophisticated regulatory mechanism in rice, which is highly possible by regulating protein expression and protein functions. Further analysis is essential to discover more regulation events to study response to PAMPs.We constructed connections between genetic control and metabolites in glycolysis/gluconeogenesis, phenylpropanoid biosynthesis and flavonoid biosynthesis. Based on that, we not only suggested a convergent response to both PAMPs in these pathways in rice but also investigated detailed regulatory mechanism that may cause the differential responses, particularly via PTMs such as ubiquitination and acetylation. By mapping the variations in metabolites and regulatory events detected in multi‐omics analyses, we systematically investigated connections between genetic control and chemical compounds in the glycolysis/gluconeogenesis, phenylpropanoid biosynthesis, flavonoid biosynthesis and linoleic acid‐associated pathways during PTI. Phenylpropanoid biosynthesis is well known to be involved in defence response in Arabidopsis (Dixon et al., 2002; Jones and Dangl, 2006). The increase in metabolites in the phenylpropionic acid biosynthesis pathway correlated with decreased ubiquitination level of the enzymes in this pathway, indicating a coordinated phenylpropanoid regulon exists by translational regulation. We detected mRNA and PTMs variations in enzyme [2.3.1.74] that was down‐regulated in response to both chitin and flg22, which may lead to increased levels of metabolites in flavonoid biosynthesis. It is widely known that flavonoids play a key role in immunity (Iranshahi et al., 2015; Petrussa et al., 2013), and many studies have described the antifungal potential of flavonoids (Chen et al., 2019; Iranshahi et al., 2015; Piasecka et al., 2015; Sekiya et al., 2021). Yet, it is still largely unknown how the flavonoids biosynthesis pathway is regulated in rice. Most of the genes, enzymes and compounds of the phenylpropanoid biosynthesis and flavonoid biosynthesis pathway changed significantly in response to chitin and flg22 (Figures 3C, 6d,e and 8). Induced accumulation of the products of the phenylpropanoid biosynthesis pathway correlated with enhanced transcription and higher abundance of the corresponding enzymes [6.2.1.12], [1.14.13.1] (Figure 8). Increases in metabolites of the phenylpropanoid biosynthesis pathway were also associated with significantly lowered ubiquitination levels of the enzymes, indicating that a coordinated phenylpropanoid regulon exists by translational regulation. We detected mRNA/proteins and PTMs variations in enzyme [2.3.1.74] that is down‐regulated in response to both chitin and flg22, which may lead to increased levels of metabolites in flavonoid biosynthesis.We highlighted linoleic acid‐associated pathways probably as an flg22‐specific responsive pathway by multiple‐omics analysis and uncovered a number of genes that may negatively modulate these pathways. As a specific response to flg22, we found linoleic acid metabolism was probably down‐regulated as almost all the detected metabolites were decreased. Linoleic acid metabolism and α‐linolenic acid metabolism pathways were specifically enriched among altered compounds and DEGs in flg22‐treated samples. As a result, induction of the linoleic acid metabolism and α‐linolenic acid metabolism will lead to accumulation of linoleic acid and α‐linolenic acid, which have been reported to be involved in callose accumulation, a classic plant PTI response (Nehela and Killiny, 2020; Zhao et al., 2015). The linoleic acid and α‐linolenic acid may be related to activation of JA‐mediated signalling pathway (Kachroo and Kachroo, 2009). It is interesting to reveal detailed roles and molecular regulatory mechanisms of linoleic acid metabolism and α‐linolenic acid metabolism, as well as linoleic acid and α‐linolenic acid, in rice PTI defence response. Moreover, we found some other pathways related to lipids (non‐saturated fatty acids biosynthesis), terpenoids (diterpenoid, sesquiterpenoid and triterpenoid biosynthesis) and vitamins (biotin metabolism and riboflavin metabolism) seem to be also specifically responsive to flg22, which would be interesting for further experimental study in the future. In fact, complex lipid changes have been widely found to produce potential molecular signatures to regulate transcriptional networks during PTI (Lim et al., 2017; Walley et al., 2013; Kachroo and Kachroo, 2009). Future studies need to focus on elucidating the functional significance of the observed specific metabolite changes and integrate all the functional data into predictive models towards holistic understanding of the signalling and metabolic processes underlying plant immunity.
Methods
Plant materials and treatment of PAMPs
Rice (Nipponbare) seedlings were prepared by cultivating the rice seeds in half strength Murashige and Skoog medium for 1 week at 28°C. Then, the rice seedlings were treated by immersing them in 8 nm chitin (Macklin, China) or 1 μm flg22 (Santa Cruz, CA) solution for different times (Chen et al., 2018). The rice seedlings treated with water were used as control. For immunoblotting experiment (0, 0.5, 1, 2, 3 and 6 h) in Figure 1B and gene expression experiments (3 h) in Figure 1C, three biological replicates were used for assessment. For transcriptome, proteome, ubiquitome and acetylome experiments, the rice seedlings were harvested after treatment of water, 8 nm chitin or 1 μm flg22 for 3 h. For each treatment, five plants were used in each of the three biological replicates (3 × 5 = 15 rice plants/treatment), and the treated rice seedlings were harvested at 3 h post‐treatment (hpt) using liquid nitrogen and stored at −70°C for further analysis. For metabolome analysis, five plants in each of the ten biological replicates (10 × 5 = 50 rice plants/treatment) were used.
Quantitative reverse transcription PCR (qRT‐PCR)
Total RNA was extracted by using TRIzol reagent (Invitrogen) from rice seedlings treated with water, chitin or flg22 for 3 h. The total RNA was treated with DNaseI (Invitrogen) to remove genomic DNA and, then, was used to synthesize the cDNA using the Reverse Transcription System (Promega). The qRT‐PCR reaction system was prepared using SYBR Green mix (TAKARA, Dalian, China), and was performed on an ImyiQ2 real‐time PCR detection system (Bio‐Rad). Expression of genes was normalized by comparing with a ubiquitin‐encoding gene OsUBQ1. Expression of the defence‐related genes in response to chitin or flg22 treatment was compared to that of water treatment for 3 h.
Detection of ROS accumulation
To detect ROS accumulation, 4‐mm‐length leaf segments were cut from 1‐week‐old rice seedlings and incubated in distilled water for 12 h. Three leaf fragments were treated with either 8 nm chitin or 1 μm flg22 in 100 mL of luminol solution, Immun‐Star HP substrate (Bio‐Rad) and 1 mL of horseradish peroxidase streptavidin (Jackson Immunoresearch). Luminescence was recorded every 20 s on a Glomax 20/20 Luminometer (Promega) for 20 min.
Western blot analysis
To detect phosphorylation level of MAPK3/6, total proteins were extracted from different samples. Rice seedling samples were grinded in liquid nitrogen, and was added to extraction buffer (50 mm HEPES (pH 7.4), 5 mm EDTA, 5 mm EGTA, 50 mm b‐glycerophosphate, 10 mm Na3VO4, 10 mm NaF, 2 mm DTT, protease inhibitor cocktail, 1% SDS or Triton X‐100). Around 50 μg total proteins were separated by SDS‐PAGE and blotted onto the PVDF membrane. The Anti‐Phospho‐p44/42 antibody (1:1000, Cell Signaling Technology) was used for immunoblotting.
Metabolites profiling
Frozen rice seedling samples were suspended with 1 mL ethanol: acetonitrile: water (2:2:1, v/v) per 80 mg of tissue fresh weight. The samples were then sonicated for 30 min, and then centrifuged for 15 min (13 000 rpm, 4°C). The aliquots were kept at −80°C and used to be analysed by ultra‐high‐performance liquid chromatography equipped with quadrupole time‐off light mass spectrometry (UPLC‐Q‐TOF/MS). Untargeted metabolic profiling was performed on an Agilent 1290 Infinity LC system (Agilent Technologies, Santa‐Clara, California). Chromatographic separation was implemented on ACQUITY HSS T3 1.8 µm (2.1 × 100 mm) columns. Detection condition was set at 25°C, with 0.3 mL/min flow rate and 2 μL sample size per each run. The two mobile‐phase solvents were water, 25 mm ammonium acetate and 25 mm ammonia (Solvent A) and acetonitrile (Solvent B). To optimize compound detection, the gradient elution programme ran as following: 95% Solvent B for 0–1 min, and linearly changed Solvent B from 95% to 65% for 1–14 min, from 65% to 40% for 14–16 min, maintained at 40% for 16–18min, linearly changed 40% to 95% for 18–18.1 min and maintained at 95% for 18.1–23 min. The entire operation was running under 4°C. To ensure the reproducibility and stability of the analysis, quality control samples were used by pooling 10 μL of each sample and analysed along with other samples. The QC samples were introduced in every 10 samples.TOF/MS was performed on positive ion mode and negative ion mode. Electrospray ionization (ESI) source conditions on triple TOF were set as following: ion source gas 2, 60 psi; curtain gas, 30 psi; source temperature, 600 °C and ionspray voltage floating, 5500 V (+) and −5500 V (−). The TOF/MS scan m/z range is 60–1000 Da, the product ion scan m/z range is 25–1000 Da, the TOF MS scan accumulation time is 0.20 s/spectra, and the product ion scan accumulation time is 0.05 s/spectra. Information‐dependent acquisition (IDA), an artificial intelligence‐based product ion scan mode, were used to detect and identify MS/MS spectra. The parameters were set as follows: declustering potential, 50 V (+) and −50 V (−); collision energy, 50 V (+) and −20 V (−); exclude isotopes within 4 Da and candidate ions to monitor per cycle: 6. The raw UPLC‐TOF/MS output files were processed to acquisite metabolite peaks by estimating signal tensity by the XCMS‐CAMERA mass scan data pipeline. The features with more than 1000 m/z were filtered. Only the peaks detected in >90% samples were kept for further analysis. We then incorporated raw data with XC‐MS workflow (Smith et al., 2006) and produced analysis by discarding features with missing values as intensity in more than 10% of total samples, after normalization and quality control. MetaAnalyst 4.0 was used to perform Pearson correlation coefficient between the replicates, and normalization by auto‐scaled method. A cluster analysis and pair‐wise comparison using MetaAnalyst were performed to determine altered metabolites between chitin‐ and flg22‐treated samples and ck samples.
Transcriptome profiling
Total RNA from rice seedling samples was extracted using Trizol regent. After quality control using Bioanalyser Agilent 2100, mRNA was enriched by PloyA and sequenced using the Illumina HiSeq‐2500 (Novogene, Beijing). Low‐quality reads and adaptor sequences were filtered using Trimmomatic 0.36. Transcript abundance was generated by Kallisto using reference genome (Bray et al., 2016). The transcript and gene read count were optimized from TPM counts using gene length variations across samples by tximport (Soneson et al., 2015). Sequencing depth, library sizes, reads dispersion and dropout were estimated by estimateParam in PowsimR package (Nedergaard et al., 2018). DESeq dataset from tximport function was used to build DESeq set for further analysis. To assess repeatability of three replicates, principal component analysis and Euclidean distance analysis were performed. To estimate and remove batch effects between replicates in each group, the sva package was used (Leek et al., 2012). Normalized transcript abundances were estimated and differentially expressed genes analysis was performed by DESeq function in DESeq2 package. Only the genes showing differential expression with log2|FC| > 1 with <0.05 P‐adjust values were defined as differentially expressed genes in pair‐wise comparisons. To perform enrichment analysis of pathways in rice, function Enrich KEGG in the R package clusterProfilter v3.16 was used (Yu et al., 2012) by inputting list of genes for corresponding analysis. For gene ontology enrichment analysis, AgriGO toolkits (Du et al., 2010) were used with the Oryza sativa Japonica database. The interested GO terms were kept using P‐adjust values cut‐offs (0.05) for significances. To visualize the network of gene ontology terms, Cytoscape 3.8 software with ClueGO plugin was used. For mapping the metabolites and transcript into pathways, Pathview (Version 3.11) (Luo et al., 2017) was used to produce detailed mapping for the selected pathways. All the figures and heatmaps were produced by ggplot2 and pheatmap package in R, unless otherwise specified.
iTRAQ proteomics profiling
Total proteins of the rice seedling samples were extracted as described by Chen et al., 2018. The samples were frozen in liquid nitrogen and ground with a pestle and mortar. Five times volume of TCA/acetone (1:9) was added to the powder and mixed by vortex. The mixture was placed at −20 °C for 4 h and centrifuged at 6000 g for 40 min at 4 °. The supernatant was discarded. Pre‐cooled acetone was added and washed for three times. The precipitation was air dried. Thirty times volume of SDT buffer was added to 30 mg powder, mixed and boiled for 5 min. The lysate was sonicated and then boiled for 15 min. After centrifuged at 14 000 g for 40 min, the supernatant was filtered with 0.22 µm filters. The filtrate was quantified with the BCA Protein Assay Kit (Bio Rad). The sample was stored at −80 °C; 20 µg of proteins for each sample, 20 µg of proteins was mixed with 5× loading buffer, and boiled for 5 min. The proteins were separated on 12.5% SDS PAGE gel (constant current 14 mA, 90 min). Protein bands were visualized by Coomassie Blue R 250 staining.A total of 200 μg proteins for each sample were dissolved into 30 μL SDT buffer (4% SDS, 100 mm DTT, 150 mm Tris HCl pH 8.0). The detergent, DTT and other low molecular weight components were removed using UA buffer (8 m Urea, 150 mm Tris HCl pH 8.0) by repeated ultrafiltration (Microcon units, 10 kD). Then, 100 μL iodoacetamide (100 mm) IAA in UA buffer was added to block reduced cysteine residues and the samples were incubated for 30 min in darkness. The filters were washed with 100 μL UA buffer three times and, then, 100 μL D isolation buffer (DS twice). Finally, the protein suspensions were digested with 4 μg trypsin (Promega) in 40 μL DS buffer overnight at 37 °C, and the resulting peptides were collected as a filtrate. The peptides of each sample were desalted on C18 Cartridges (Empore™ SPE Cartridges C18, standard density), bed I.D. 7 mm, volume 3 mL; Sigma‐Aldrich, St. Louis, MO), concentrated by vacuum centrifugation and reconstituted in 40 µL of 0.1% (v/v) formic acid. The peptide content was estimated by UV light spectral density at 280 nm using an extinctions coefficient of 1.1 of 0.1% (g/L) solution that was calculated on the basis of the frequency of tryptophan and tyrosine in vertebrate proteins.A total of 100 μg peptide mixture of each sample was labelled using iTRAQ reagent according to the manufacturer’s instructions (Applied Biosystems, Foster City, CA). iTRAQ labelled peptides were fractionated by SCX chromatography using the AKTA Purifier System (GE Healthcare, Boston, MA). The dried peptide mixture was reconstituted and acidified with buffer A (10 mm KH2PO4 in 25% of ACN, pH 3.0) and loaded onto a PolySULFOETHYL 4.6 × 100 mm column (5 μm, 200 Å, PolyLC Inc., Columbia, MD). The peptides were eluted at a flow rate of 1 mL/min with a gradient of 0%–8% buffer B (500 mm KCl, 10 mm KH2PO4 in 25% of ACN, pH 3.0) for 22 min, 8–52% buffer B during 22–47 min, 52%–100% buffer B during 47–50 min, 100% buffer B during 50–58 min and buffer B was reset to 0% after 58 min. The elution was monitored by absorbance at 214 nm, and fractions were collected every 1 min. The collected fractions were desalted on C18 Cartridges (Empore™ SPE Cartridges C18, standard density), bed I.D. 7 mm, volume 3 mL (Sigma) and concentrated by vacuum centrifugation.Each fraction was injected for nano‐LC‐MS/MS analysis. The peptide mixture was loaded onto a reverse phase trap column (Thermo Scientific Acclaim PepMap100), 100 μm × 2 cm and nanoViper C18 connected to the C18 reversed phase analytical column (Thermo Scientific Easy Column, 10 cm long, 75 μm inner diameter and 3 μm resin) in buffer A (0.1% Formic acid) and separated with a linear gradient of buffer B (84% acetonitrile and 0.1% formic acid) at a flow rate of 300 nL/min controlled by IntelliFlow technology.LC‐MS/MS analysis was performed on a Q Exactive mass spectrometer (Thermo Scientific, Waltham, MA) that was coupled to Easy nLC (Proxeon Biosystems, now Thermo Fisher Scientific) for 120 min. The mass spectrometer was operated in positive ion mode. MS data were acquired using a data‐dependent method dynamically choosing the top 10 most abundant precursor ions from the survey scan (300–1800 m/z) for HCD fragmentation. Automatic gain control (AGC) target was set to 3e6, and maximum inject time to 10 ms. Dynamic exclusion duration was 40.0 s. Survey scans were acquired at a resolution of 70 000 at m/z 200 and resolution for HCD spectra was set to 17,500 at m/z 200 and isolation width was 2 m/z. Normalized collision energy was 30 eV and the underfill ratio, which specifies the minimum percentage of the target value likely to be reached at maximum fill time, was defined as 0.1%. The instrument was run with peptide recognition mode enabled.MS/MS spectra were searched using the MASCOT engine (Matrix Science, London, UK; version 2.2) embedded into Proteome Discoverer 1.4 using 20 ppm Mass Torrance. The missing values in the dataset were imputed use DEP package in R, and the peptide ratios were normalized by the median protein ratio. Correlation analysis and principal component analysis were performed using R. For differentially expressed proteins analysis, t‐test was performed by each pair‐wise comparison and P‐values were corrected by Bonferroni correction. Differentially expressed proteins were defined by log2‐fold changes ≥1, with <0.05 P‐adjust values in each pair‐wise comparison.
The proteins analysis for PTMs
The dataset for ubiquitination was used from previous study (Chen et al., 2018). To identify peptides modified by acetylation, 10 mg of each group of samples were used to add 45 mL TCA/acetone (1:9, containing 65 mm DTT). The samples were precipitated overnight at −20 °C and then sent for centrifuge at 7000 ×g for 20 min. The supernatant was discarded, following the addition of 40 mL of pre‐cooled acetone with 7000 ×g centrifuge for 15 min to remove the supernatant. The washing process was repeated three times and the precipitate was dried at room temperature. Take the dried powder, add an appropriate amount of SDT lysis solution, boil in a water bath for 15 min and repeat the boiling twice. Ultrasound the sample (100w, ultrasonic 10s, interval 10s and repeat 10 times). Centrifuge at 13 400 rpm for 20 min and collect the supernatant. Bicinchoninic Acid (BCA) method was used for protein concentration determination. After that, 20 μg of each protein sample was added with 5× loading buffer and was incubated in boiling water bath for 5 min. The samples were centrifuged at 14 000 ×g for 10 min to take the supernatant and performed 12.5% SDS PAGE electrophoresis. For enzymatic hydrolysis of protein, 10 mg of each sample was added with DTT to the final concentration of 10 mm at 37 °C for 2.5 h, and then cooled to the room temperature. Add IAA to a final concentration of 50 mm and avoid from light for 30 min. Each sample was added five times the volume of water to dilute the UA concentration to 1.5 m, following addition of trypsin at a ratio of 1:50 and digest for 18h at 37 °C. For enrichment of acetylated peptides, each sample was lyophilized by adding 1.4 mL of pre‐chilled IAP Buffer, then adding pre‐treated anti‐Ac‐K antibody beads (PTMScan, Acetyl Lysine Motif Ac K Kit, Cell Signaling Technology, Beverly, MA). The samples were incubated at 4 °C for 1.5 h, then centrifuged at 2000 g for 30 s, and the supernatant was discarded. Anti‐Ac‐K antibody beads are washed three times with 1 mL of pre‐cooled IAP Buffer, and then washed three times with 1 mL of pre‐cooled water. The clear samples were added with 40 μL of 0.15% TFA to the washed anti‐Ac‐K antibody beads at room temperature for 10 min, following addition of a total of 0.15% TFA twice. Centrifuge at 2000 g for 30S, take the supernatant and desalinate by C18 STAGE Tips.The samples were processed for the enriched acetylated peptides for LC‐MS/MS analysis. The nano‐litre flow rate HPLC liquid system, System EASY nLC1000, was used for separation. Liquid A was 0.1% formic acid acetonitrile aqueous solution and 2% acetonitrile; and liquid B was 0.1% formic acid acetonitrile aqueous solution and acetonitrile is 84%. Chromatographic column Thermo EASY column SC200 7/11 150 μm*100 mm was balanced to100% AA liquid. The samples were loaded by the autosampler to liquid equilibrium, onto the Thermo EASY column SC001 traps 150 μm*20 mm (RPEASY column SC001 traps 150 μm*20 mm (RP‐CC1818)) ((ThermoThermo)) by the autosampler, and then passed through the chromatographic column and separated by the chromatographic column. The relevant liquid‐phase gradients are as follows: For 0 to 110 min, B liquid was increased linearly from 0% to 55%; for 110 to 118 min, Bliquid was increased linearly from 55% to 100% and for 118 to 120 min, B liquid was maintained at 100%. Capillary high‐performance liquid chromatography was used after separation. The Q Exactive mass spectrometer (Thermo Finnigan) was used for mass spectrometry. Analysis time: 120 min, detection method: positive ion and parent ion scan range: 350–1800 m/z. For peptides and peptides, the mass‐to‐charge ratio of the fragments was collected according to the following method: 20 pieces are collected after each full scan, Fragment Map MS 2 scan HCD. MS 1 has a resolution of 70 000 at M/Z 200 MS 2 at M/Z. The resolution was 17 500 at 200 h. The raw data were imported into Maxquant (Version 1.3.0.5), and annotated peptides using amino acid FASTA file obtained from NCBI, by using parameters as below: main search ppm: 6; missed cleavage: 4; MS/MS tolerance pmm: 20; peptide FDR:0.01; protein FDR:0.01 and fixed modification carbamidomethyl (C).
Conflict of interest
The authors declare no conflicts of interest.
Author contributions
W.L. and X.‐L.C. designed the project. C.L., Z.L., and X.Z. performed the experiments. B.T., X.‐L.C., W.L., and G.‐L.W. performed the multi‐omics data analysis and interpretation of the results. S.Z. assisted in data analysis and manuscript preparation. W.L., X.‐L.C. and B.T. wrote the manuscript.Figure S1 Overview of metabolomics analysisFigure S2 Overview of transcriptomics analysisFigure S3 Mapping of cellular processes and metabolic pathways induced by chitin and flg22 treatmentFigure S4 Overview of proteomics analysisFigure S5 Overview of ubiquitomics analysisFigure S6 Overview of acetylomics analysisFigure S7 Integrative analysis of multi‐omics data of rice affected by chitin‐ and flg22‐treatmentFigure S8 Mapping of cellular processes and metabolic pathways suppressed by chitin and flg22 treatmentClick here for additional data file.Dateset S1 All features from metabolomics.Click here for additional data file.Dateset S2 Significantly altered chemical compounds upon chitin and flg22 induction.Click here for additional data file.Dateset S3 Metabolic pathways enriched by altered compounds detected in metabolomics.Click here for additional data file.Dateset S4 Differentially expressed genes in rice responding to chitin and flg22.Click here for additional data file.Dateset S5 Differentially expressed genes related to rice innate immunity.Click here for additional data file.Dateset S6 Expression of all proteins detected by iTRAQ analysis.Click here for additional data file.Dateset S7 Differentially expressed proteins in rice responding to chitin and flg22.Click here for additional data file.Dateset S8 Expression of ubiquitinated proteins detected in ubiquitome.Click here for additional data file.Dateset S9 Quantitative analysis of differentially ubiquitinated proteins.Click here for additional data file.Dateset S10 Expression of acetylated proteins detected in acetylome.Click here for additional data file.Dateset S11 Quantitative analysis of differentially acetylated proteins.Click here for additional data file.
Authors: Damian Szklarczyk; Alberto Santos; Christian von Mering; Lars Juhl Jensen; Peer Bork; Michael Kuhn Journal: Nucleic Acids Res Date: 2015-11-20 Impact factor: 16.971