Scientific Papers

Transcriptome and metabolome analyses revealed the response mechanism of pepper roots to Phytophthora capsici infection | BMC Genomics

Phenotypic differences between A198 and A204 in response to P. capsici infection

The two pepper cultivars used in this study were obtained by screening a core collection of 200 pepper germplasm resources for P. capsici resistance (Fig. 1a). Based on disease phenotypes, this study selected two genotypes with markedly different resistance levels. A204, the resistant genotype, showed the same resistance as the resistant ‘Criollo de Morelos 334’ (CM334) which is one of the most promising sources of resistance to P. capsici in pepper [40], and had no disease symptoms at 96 hpi with P. capsici (Fig. 1b). For susceptible genotype A198, all plants were wilted and had a 2–3 cm constriction at the stem base, with defoliation of leaves at 96 hpi (Fig. 1b). To determine the appropriate sampling time points for RNA sequencing (RNA-seq) and metabolomic studies, continuous observation of the symptoms of the susceptible genotype post-inoculation was performed. There were no apparent symptoms at 12 hpi, while brown lesions formed on the junction of the root and stem at 36 hpi, and about 0.5 and 1 cm constrictions were present at the stem base at 60 and 72 hpi, respectively (Fig. 1c). Therefore, 24 hpi was considered a suitable time point for the early stage of infection, which was to be in the biotrophic phase, and 48 hpi, the mid-late stage, was to be in the necrotrophic phase.

Fig. 1
figure 1

Phenotypic observation of pepper inoculated with P. capsici. (a) The phenotype of partial pepper germplasm resources 10 days after inoculation for the screening of Phytophthora root rot resistance. (b) The phenotype of the susceptible genotype A198 and the resistant genotypes A204 and CM334 at 96 h post-inoculation (hpi). (c) Symptoms of susceptible genotype A198 at 12, 36, 60, and 72 hpi. OV means the overall view. PV represents the partial view

Overview of pepper root transcriptomic responses to P. capsici infection

Transcriptome analysis was performed on the pepper roots of the two pepper genotypes at 0, 24, and 48 hpi. A total of 118.08 Gb clean data were obtained after removing reads containing adapters, ploy-N, and low-quality reads. The clean data obtained for each sample reached 5.80 Gb, and the Q20 and Q30 base percentage were at least 97.91 and 93.97%, respectively (Table 1). The alignment results showed that at least 88.69% of the high-quality reads were successfully mapped to the reference genome of CM334. Of the mapped reads, 84.71–88.33% were mapped to exon regions, 4.14–5.72% to intron regions, 7.39–9.39% to intergenic regions, and 0.14–0.18% to spliced regions (Supplementary Materials 1:Table S1), which implies that a set number of transcripts were possibly derived from alternative mRNA splicing or new genes. The Pearson’s correlation coefficients (R2) between biological replicates were all above 0.99, demonstrating a high degree of biological reproducibility among the samples (Supplementary Materials 7: Fig. S1a). Principal component analysis (PCA) of all transcripts showed that the first two principal components (PCs) explained 74.8% of the total variation, and differences between the genotypes could be clearly displayed, regardless of whether they were infected or not (Supplementary Materials 7: Fig. S1b). Upon P. capsici infection, transcriptional profiles changed in opposite directions between the two genotypes, indicating that the two genotypes may resist P. capsici infection through different response mechanisms.

Table 1 General information of sequencing reads and reads that mapped to the reference genome

Functional annotation and identification of DEGs

Through database search of detected transcripts, 35,642 genes (including 4750 new genes) were annotated, of which 9546, 25,863, 23,776, 17,636, 25,893, 21,600, 29,437, and 35,576 genes were annotated in the Cluster of Orthologous Group (COG), Gene Ontology (GO), KEGG, Karyotic Ortholog Groups (KOG), Pfam, Swiss-Prot, evolutionary genealogy of genes: Non-supervised Orthologous Groups (eggNOG), and non-redundant (NR) databases, respectively. To identify DEGs, pairwise comparisons were independently performed between each time point during infection and at 0 hpi in both genotypes. A total of 1517 and 1227 DEGs were identified in the resistant genotype, and 1057 and 1968 DEGs were determined in susceptible genotype A198 at 24 and 48 hpi, respectively (Fig. 2; Supplementary Materials 2: Table S2). The DEGs from the same genotype overlapped more, but the DEGs from different genotypes were mostly unique (Fig. 2a). Among the DEGs, the number of upregulated DEGs was higher than that of downregulated DEGs, and the downregulated DEGs in A204 were far more numerous than those in A198 at all time points (Fig. 2b). Additionally, a hierarchical cluster analysis of all DEGs was performed to assess the reproducibility of RNA-sequencing data, and the cluster heatmap showed that the three biological replicates of each sample group were clustered into one cluster, which indicated the reliable reproducibility of RNA-sequencing (Supplementary Materials 7: Fig. S2).

Fig. 2
figure 2

Transcriptional changes upon P. capsici infection in A198 and A204. (a) Venn diagram of differentially expressed genes (DEGs). R1/ R2 vs. R0 represent a comparison between 24/ 48 h post-inoculation (hpi) and 0 hpi, in the resistant genotype. S1/ S2 vs. S0 represent a comparison between 24/ 48 hpi and 0 hpi in the susceptible genotype. (b) Number of up- and downregulated DEGs of two genotypes at 24 and 48 hpi

DEGs enrichment analysis

To determine the pathways of the DEGs involved in pepper response to P. capsici infection, a KEGG enrichment analysis of the DEGs was performed. A total of 37 and 26 KEGG pathway categories were enriched in the up- and downregulated DEGs, respectively (Fig. 3). Among them, the 13 enriched KEGG pathways shared by both up- and downregulated DEGs mainly involved carbohydrate metabolism, energy metabolism, signal transduction, and the biosynthesis of other secondary metabolites. Correspondingly, 24 and 12 KEGG pathway categories were uniquely enriched in up- and downregulated DEGs, respectively (Fig. 3). For signal transduction, the “plant hormone signal transduction pathway” and “MAPK signaling pathway” were significantly enriched in the susceptible genotype at 24 and 48 hpi among the upregulated DEGs. The number of DEGs involved in the two pathway categories showed an increasing trend from 24 to 48 hpi in the susceptible genotype, while the opposite trend was found in the resistant genotype. These results suggest that the resistant genotype may respond to P. capsici earlier than the susceptible genotype. Notably, most pathways involved in the “metabolism of terpenoids and polyketides” and “biosynthesis of other secondary metabolites” were significantly enriched in the resistant genotype at 24 hpi, including “sesquiterpenoid and triterpenoid biosynthesis”, “diterpenoid biosynthesis”, “carotenoid biosynthesis”, “phenylpropanoid biosynthesis”, “flavonoid biosynthesis”, “anthocyanin biosynthesis”, and “flavonoid and flavonol biosynthesis”, which were considered to be related to the response to biotic and abiotic stress. In terms energy metabolism, namely “photosynthesis-antenna proteins”, “carbon fixation in photosynthetic organisms”, and “nitrogen metabolism” were significantly enriched in the susceptible genotype in the upregulated DEG group, while these pathways were significantly enriched in the resistant genotype in the downregulated DEG group. Unexpectedly, most pathways involved in amino acid metabolism, such as “tyrosine metabolism”, “cysteine and methionine metabolism”, “valine, leucine and isoleucine degradation”, “arginine and proline metabolism”, and “phenylalanine metabolism”, were significantly enriched only in the susceptible genotype. The same result was found for the pathways related to “lipid metabolism” and “fatty degradation”. These results may suggest that the enhanced energy metabolism, as well as amino acid and lipid catabolism, was enhanced to meet energy requirements in the susceptible genotype.

Fig. 3
figure 3

Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis of up- (a) and downregulated (b) differentially expressed genes (DEGs). The number in the box is the number of DEGs enriched in that pathway per sample

DEGs associated with disease resistance

To better understand the network of DEGs that responded to P. capsici infection, the transcriptional changes of resistant and susceptible pepper roots were visualized using MapMan software. The heatmap shows a general overview of the differences in the response to P. capsici infection between the two genotypes (Supplementary Materials 7: Fig. S3). As already highlighted by the DEG analysis, there were more DEGs involved in cell wall strengthening, signaling, and secondary metabolite synthesis at 24 hpi in the resistant genotype than in the susceptible genotype (Supplementary Materials 7: Fig. S3a and b), while at 48 hpi, the opposite response to P. capsici infection was observed in the DEGs of the two genotypes (Supplementary Materials 7: Fig. S3c and d), indicating that the resistant genotype reacted more quickly to infection than the susceptible genotype.

Structural defense. The plant cell wall is the first line of defense, providing the plant with initial defense and signal perception against pathogen attack. In this study, five subtilisin-like protease genes were found to be especially induced in the resistant genotype at 24 hpi (Supplementary Materials 3: Table S3). Interestingly, from 24 to 48 hpi, the number of upregulated DEGs encoding for the LRR protein in the resistant genotype decreased, while the number of DEGs in the susceptible genotype continued to increase. Furthermore, genes responsible for cell wall strengthening, such as proline-rich protein (PRP), cellulose synthase (CesA), and syntaxin, were also differentially expressed. Unexpectedly, the strong repression of most unigenes encoding proline-rich protein occurred in both genotypes at 24 and 48 hpi, while the comparison of the baseline levels of the detected DEGs between A204 and A198 (R0 vs. S0) showed higher levels of six proline-rich protein gene transcripts (Supplementary Materials 3: Table S3), which may contribute to improved pathogen defense in the resistant genotype. For cellulose synthase, two genes were upregulated in the resistant plants at 48 hpi. Only one gene, LOC107853789, was upregulated in the susceptible plants, and this gene had a higher baseline level (25.35-fold) in the resistant genotype (Supplementary Materials 3: Table S3). In addition, the A204 response to P. capsici was characterized by the earlier induction of cell wall degradation-related enzymes, such as polygalacturonase (PG) and pectinesterase (PE), as well as pectate lyase. At 24 hpi, six DEGs encoding PG were induced in A204, but only one was induced in the susceptible genotype. More DEGs coding for PE and pectate lyase were upregulated at 48 hpi in A198. These results indicate that the degradation of the cell wall, as well as the response to pathogen attack, occurred earlier in the resistant genotype than in the susceptible genotype.

Signal transduction. To mount an effective defense, plants rapidly transmit stress signals to trigger various downstream defense mechanisms through different signaling pathways, including calmodulin (CaM) and phytohormone signaling pathways. Nine CaM-related DEGs were upregulated in the resistant genotype (Supplementary Materials 3: Table S3), suggesting a key role of the Ca2+-dependent signaling pathway in resistant genotype A204. It is well-known that the SA, JA, and ET pathways play important roles in signal transduction in response to biotic stress. In this regard, three SA-related DEGs encoding ankyrin repeat-containing proteins were induced at 24 hpi in A204. A total of 37 DEGs related to ET-responsive factors were identified. among them, 18 DEGs were expressed in A198 and nine in A204 (Supplementary Materials 3: Table S3). For JA-related genes, three DEGs coding for lipoxygenase (LOX), one for a ZIM domain-containing protein, and one for TIFY were upregulated in the susceptible genotype, but only one DEG encoding TIFY was induced in A204. Furthermore, 11 auxin-related genes were identified in the two genotypes. Interestingly, they were differentially expressed in A198 only at 48 hpi but differentially expressed at 24 and 48 hpi in A204. Similar to the changes in auxin-related DEGs, three abscisic acid (ABA) receptor genes were upregulated in A198 only at 48 hpi. In addition, activation of the MAPK cascade and WRKY transcription factor family members were also observed (Supplementary Materials 7: Fig. S3). MPKK5 genes were suppressed in both genotypes at 48 hpi. Six and seven WRKY genes were induced in A198 at 24 and 48 hpi, respectively (Supplementary Materials 3: Table S3). These results suggest that the two genotypes may transmit stress signals through different signaling pathways.

Chemical defense. When a pathogen breaks through the first line of defense, plants usually fight pathogens and repair themselves through the employment of numerous specialized proteins, such as xylem proteinase and pathogenesis-related protein (PR). As shown in Fig. S3 (Supplementary Materials 7), more genes encoding specialized proteins were differentially induced in A204 at 24 hpi, including xylem cysteine proteinase 1 (XCP1), pathogenesis-related protein 1 (PR1), endo-1,3(4)-beta-glucanase (PR2), chitinase (PR3 and PR4), thaumatin-like protein (PR5), and peroxidase (PR9) (Supplementary Materials 3: Table S3). Two XCP1 genes (LOC107851745 and LOC107866767) were uniquely induced in A204. For PR1, two DEGs were upregulated only at 24 hpi in the resistant genotype, while one was upregulated in the susceptible genotype at 48 hpi, suggesting a late defensive response in the susceptible materials. Three PR2-related genes, namely LOC107866197, LOC107839367, and newGene_10657, were upregulated only in A204, while one (newGene_10487) was induced in A198 (Supplementary Materials 3: Table S3). The PR3-, PR4-, and PR5-related DEGs were most induced in the susceptible genotype, especially nine chitinase-related genes that were upregulated at 48 hpi (Supplementary Materials 3: Table S3), which demonstrated that a strong defensive response occurred at that time. In addition, eight genes implicated in biosynthesis of flavonoid were uniquely upregulated in A204 at 24 hpi (Supplementary Materials 3: Table S3).

Validation of the expression of selected DEGs

To confirm the validity of the transcriptome data, real-time quantitative polymerase chain reaction (RT-qPCR) analysis was performed on 12 selected genes (Fig. 4). These selected genes were involved in signal transduction (CaM LOC107868363, CaMBP LOC107862193, and SBT LOC107872604), cell wall enhancement (PRP LOC107852566 and CesA LOC107853789), pathogen resistance (PR2 LOC107866197, PR3 LOC107859802, and PR4 LOC107840165), and flavonoid biosynthesis (CHS LOC107850995, CHI LOC107852750, F3H LOC107859880, and FLS LOC107876027), which were either directly or indirectly linked to plant resistance. The comparison between the two techniques revealed substantial agreement for all 12 genes differentially expressed upon pathogen infection.

Fig. 4
figure 4

Comparison of real-time quantitative polymerase chain reaction (RT-qPCR) and RNA-Sequencing results. R1/ R2 vs. R0 represent a comparison between 24/ 48 h post-inoculation (hpi) and 0 hpi in the resistant genotype. S1/ S2 vs. S0 represent a comparison between 24/ 48 hpi and 0 hpi in the susceptible genotype

Metabolite profiles

Untargeted metabolomics analyses were performed on pepper root extracts using an UPLC-MS/MS platform. After data filtering and identification, 688 compounds (including 514 in positive ion mode and 174 in negative ion mode) were acquired from the root extracts of all samples (Supplementary Materials 4: Table S4). PCA of the metabolome data showed tight clustering of the replicate samples of both genotypes and quality control samples, confirming the reproducibility of the results (Fig. 5a). PCA score plots further revealed that, similar to the transcriptome data, differences between genotypes could be clearly displayed based on PC1, regardless of whether they were infected. A heatmap of all analyzed ions revealed distinct hierarchical clustering of the samples based on genotype, and ions differentially accumulated with either different genotypes or different time points (Fig. 5b). To identify the DAMs, orthogonal partial least squares-discriminant analysis (OPLS-DA) models were constructed to maximize the difference between the experimental sample groups (R1, R2/S1, and S2) and the control sample groups (R0/S0). The OPLS-DA models (Q2 > 0.9) showed reliable predictability and significant biochemical perturbation in the experimental groups (Supplementary Materials 7: Fig. S4 and 5).

Metabolites with fold-changes of ≥ 2.0 or ≤ 0.5, VIP values > 1.0, and p-values < 0.05 were selected as DAMs. There were 53 and 80 DAMs in the resistant genotype at 24 and 48 hpi, respectively (Supplementary Materials 5: Table S5). The top 10 significantly altered metabolites at 24 hpi are showed in Fig. 5c. Among them, the upregulated DAMs mainly consisted of phenolic compounds, including flavonoids (neohesperidin and formononetin 7-O-glucoside) and phenylpropanoids (isoscopoletin, umbelliferone), and the top 10 downregulated DAMs mainly included lipids (sphinganine, 12,13-DiHOME, and 16-hydroxypalmitic acid). For the susceptible genotype, 32 and 108 DAMs were identified at 24 and 48 hpi, respectively (Supplementary Materials 5: Table S5). The top 10 upregulated DAMs mainly included carbohydrates (maltohexaose, cellotetraose, and raffinose) and nucleotide derivates (deoxyinosine and deoxycytidine), and the downregulated DAMs mostly consisted of nucleotide derivatives and phenylpropanoids (Fig. 5d).

Fig. 5
figure 5

Metabolic changes upon P. capsici infection in A198 and A204. (a) Principal component analysis (PCA) of all metabolites detected at 0 (R0 and S0), 24 (R1 and S1), and 48 (R2 and S2) hours post-inoculation (hpi) in both genotypes. (b) Clustering analysis heat map of the expression of differentially accumulated metabolites (DAMs) detected in each sample. Color coding represents the metabolite concentrations with high (red) and low (blue) expression. Top 10 DAMs in the resistant genotype A204 (c) and susceptible genotype A198 (d)

Association of transcriptomic and metabolomic changes involved in the response of pepper to P. capsici infection

To obtain more information on the physiological changes of pepper in response to P. capsici, this work concentrated on the connection between gene expression and metabolite changes. Based on the KEGG analysis results of DEGs and DAMs, a heatmap of the flavonoid biosynthesis pathways flow chart was drawn in this study (Fig. 6). The results showed that five flavonoid compounds, namely neohesperidin, naringin, apiin, galangin, and formononetin 7-O-glucoside, were accumulated in the resistant genotype A204, while in the susceptible genotype A198, these DAMs were decreased or unchanged. The expression patterns of most genes were similar to those of the metabolites (Fig. 6). For two chalcone synthase (CHS) genes, LOC107850995 was upregulated in A204 at 24 hpi, while the expression of LOC107871238 was suppressed in A198 at both time points (Fig. 6; Supplementary Materials 3: Table S3). Chalcone isomerase (CHI) catalyzes the cyclization of chalcone into flavanone, and the gene LOC107852750 encoding CHI was strongly induced only in A204 at 24 hpi (Fig. 6; Supplementary Materials 3: Table S3). As with CHI, five genes involved in flavonoid biosynthesis, namely flavanone 3-hydroxylase (F3H), flavonoid 3’,5’-hydroxylase (F3’5’H), dihydroflavonol reductase (DFR), flavonol-3-O-glucoside L-rhamnosyltransferase (FG2), and vestitone reductase (VR), were uniquely upregulated in resistant genotype A204 at 24 hpi. These results of combined transcriptome and metabolome analysis suggest that the resistance of genotype A204 could be related to its special ability or earlier ability to regulate gene expression in flavonoid biosynthesis pathways and the accumulation of flavonoids.

Fig. 6
figure 6

Heatmap of the log2-fold changes in genes and metabolites involved in the flavonoid biosynthesis pathways in both genotypes upon P. capsici infection at 24 and 48 h post-inoculation (hpi). In the schemes of cascades, compounds are shown in black font and genes in blue. ANS, anthocyanidin synthase; ANR, anthocyanidin reductase; CHI, chalcone isomerase; CHS, chalcone synthase; C12RT1, flavanone 7-O-glucoside 2’’-O-beta-L-rhamnosyltransferase; CYP81E, isoflavone/4’-methoxyisoflavone 2’-hydroxylase; CYP93B2, flavone synthase II; DFR, dihydroflavonol reductase; F3H, flavanone 3-hydroxylase; F3’5’H, flavonoid 3’,5’-hydroxylase; FG2, flavonol-3-O-glucoside L-rhamnosyltransferase; FLS, flavonol synthase; FNS, flavone synthase; VR, vestitone reductase. The heatmap was drawn based on the KEGG pathway maps [41]

Source link