Pathophysiology and transcriptomic responses of Pinus armandii defenses to ophiostomatoid fungi

Abstract

Pinus armandii Franch. is an ecologically and economically important evergreen tree species native to western China. Dendroctonus armandi Tsai and Li and pathogenic ophiostomatoid fungi pose substantial threats to P. armandii. With the interplay between species, the defense mechanisms of P. armandii have evolved to withstand external biotic stressors. However, the interactions between P. armandii and pathogenic ophiostomatoid fungal species/strains remain poorly understood. We aimed to analyze the pathophysiological and molecular changes in P. armandii following artificial inoculation with four ophiostomatoid species (Graphilbum parakesiyea, Leptographium qinlingense, Ophiostoma shennongense and Ophiostoma sp. 1). The study revealed that L. qinlingense produced the longest necrotic lesions, and G. parakesiyea produced the shortest. All strains induced monoterpenoid release, and monoterpene levels of P. armandii were positively correlated with fungal virulence (R² = 0.93, P < 0.01). Co-inoculation of two dominant highly (L. qinlingense) and weakly virulent (O. shennongense) pathogens reduced the pathogenicity of the highly virulent fungi. Transcriptomic analysis of P. armandii (LQ: L. qinlingense treatments, QS: co-inoculation treatments and OS: O. shennongense treatments) showed that the expression pattern of differentially expressed genes (DEGs) between QS and OS was similar, but different from that of LQ. The DEGs (LQ vs QS) involved in flavonoid biosynthesis and phenylpropanoid biosynthesis were downregulated. Notably, compared with LQ, QS significantly decreased the expression of host defense-related genes. This study provides a valuable theoretical basis for managing infestations of D. armandi and associated ophiostomatoid fungi.

Introduction

Bark beetles (Coleoptera, Curculionidae, Scolytinae) breed in conifers and form complex ectosymbiotic relationships with several fungi, especially ophiostomatoid fungi (Ophiostomataceae, Ascomycota), that can contribute to the colonization success of bark beetles (Bentz and Six 2006; Lieutier et al. 2009; DiGuistini et al. 2011; Wadke et al. 2016; Davis et al. 2019; Dzurenko and Hulcr 2022). Notably, bark beetles are commonly associated with various fungal complexes that exhibit different or overlapping ecological roles (Kirisits 2004; Zhao et al. 2019). Ophiostoma clavatum (pathogenic fungi) and Ophiostoma macrosporum (nutritional fungi) are associated with Ips acuminatus, which co-infests Scots pine (Pinus sylvestris L.) (Paine et al. 1997; Guérard et al. 2000; Harrington et al. 2010; Villari 2012); several ectosymbiotic fungi that can detoxify spruce phenolics are associated with Ips typographus (Hammerbacher et al. 2013; Zhao et al. 2019).

Fungi associated with bark beetles are typically invasive pathogenic species (Six and Wingfield 2011). The pathogenicity of fungi refers to their ability to kill hosts, and the level of pathogenicity corresponds to the virulence of the fungus. Ophiostomatoid fungi are necrotrophic pathogens, with several important tree phytopathogens substantially penetrating the xylem and phloem to form necrotic lesions or kill hosts (Kirisits 2004; Villari et al. 2012; Santini and Faccoli 2015). As with all interspecies interactions, to mitigate the predations of herbivores and pathogens, conifers have developed a multipurpose suite of biochemical defense mechanisms (i.e. oleoresin terpene and phenolic compounds) through the evolutionary arms race (Adams et al. 2011; Kushalappa et al. 2016; Celedon and Bohlmann 2019; Hammerbacher et al. 2020; de Ruiter 2022).

Terpenes stored in resin ducts in the phloem and sapwood are sophisticated defense barriers in conifers (Celedon and Bohlmann 2019; Zaman et al. 2023). The release of terpenes by trees increases significantly following a bark beetle invasion or pathogenic fungal infections (Klepzig et al. 1996; Ralph et al. 2006; Raffa 2014; Erbilgin 2019). In vitro trials have shown that conifer terpenes can directly or indirectly affect herbivore mortality as toxins (Davis and Hofstetter 2011; Zhao et al. 2011; Cale et al. 2017; Ullah et al. 2021), feeding barriers, or attractants of parasitoids and predators (Delphia et al. 2006; Raffa 2014; Fang et al. 2020). Additionally, they effectively inhibit the growth of ophiostomatoid fungi, such as Grosmannia clavigera, in a concentration-dependent manner (Cale et al. 2017; Ullah et al. 2021). Because highly virulent fungi induce longer necrotic lesions and higher monoterpene concentrations upon infection than weakly virulent fungi (Fäldt et al. 2006; Villari et al. 2012; Yamaoka 2017), the lesion length and increase in monoterpene levels have been previously used as indicators to assess the virulence of fungi following experimental inoculation.

Although phenols have received comparatively minimal attention in scientific research, they are crucial in the immune response of conifers (Hammerbacher et al. 2020). Catechins and taxifolin have inhibitory effects on beetles, and high concentrations of some flavonoids impede fungal development (Wallin 2005; Faccoli and Schlyter 2007; Hammerbacher et al. 2014, 2019).

In addition to pathophysiological responses, pathogen inoculation induces changes in the host defense mechanisms through upregulation and downregulation of gene expressions (Ganthaler et al. 2017; Nemesio-Gorriz et al. 2017). Powerful molecular tools are now available to evaluate and elucidate the mechanisms of immune system responses in conifers. The low cost and high accuracy of RNA-seq have enabled its successful application in analyzing host–pathogen interaction mechanisms and mining functional genes (Liu et al. 2017, 2023; Wilkinson et al. 2022). Recently, our knowledge of defense metabolic pathways in conifers has expanded with the application of genome sequencing and annotation of several conifer species, such as Norway spruce (Nystedt et al. 2013), Pinus tabuliformis Carrière (Niu et al. 2022) and Larix kaempferi (Lamb.) Carr (Sun et al. 2022). Defense metabolic pathways, such as terpenoid biosynthesis, flavonoid and phenylpropanoid biosynthesis and mitogen-activated protein kinase (MAPK) signaling pathways, are activated to resist invasion by pathogenic fungi (Celedon and Bohlmann 2019; Liu et al. 2022; Wilkinson et al. 2022).

In China, Pinus armandii Franch. is an evergreen tree of the Pinaceae family; it substantial contributes to economic growth and ecological development (Li et al. 2012). However, it has been substantially affected by the recent trends in global warming, droughts and the serious threat of Dendroctonus armandi Tsai and Li outbreaks, including the colonization of trees by phytopathogenic ophiostomatoid fungi associated with this bark beetle species (Cai 1980; Chen et al. 1999; Tang and Chen 1999; Zhang et al. 2015). Common ophiostomatoid fungi associated with D. armandi include Graphium parakesiyea, Graphium pseudormiticum, Leptographium qinlingense, Leptographium wushanense, Ophiostoma brevicolle, Ophiostomafloccosum, Ophiostoma shennongense and Ophiostoam sp. 1 (Wang et al. 2022). Research on the pathogenicity of ophiostomatoid fungi associated with D. armandi on P. armandii is limited. Except for studies on L. qinlingense-inoculated P. armandii (Tang and Chen 1999; Chen and Tang 2002; Pham et al. 2014), fungal pathogenicity remains largely undetermined.

In this study, we sought to analyze the pathophysiological and molecular changes in P. armandii following artificial inoculation with ophiostomatoid species. Four representative ophiostomatoid species (L. qinlingense and O. shennongense, two dominant and widespread species in the ophiostomatoid community, and G. parakesiyea and Ophiostoma sp. 1, two less frequently obtained but widely distributed ophiostomatoid fungal species across the sampling area) associated with D. armandi on P. armandii were selected for pathogenicity analysis. Chemical analyses were performed to assess changes in monoterpenoid content during host responses to pathogens to clarify the differences in virulence among different ophiostomatoid fungal species/strains. Additionally, a transcriptome assay was used to analyze differential gene expression between co-single inoculation of P. armandii with fungal species. The results were used to determine the defensive metabolic pathways of the host affected by the ophiostomatoid fungi complex responsible for reducing the pathogenicity of highly pathogenic fungi. Our findings provide insights into the pathophysiological and molecular mechanisms driving the host–fungus interaction in the P. armandii disease system, as well as providing a scientific basis for managing D. armandi and phytopathogenic ophiostomatoid fungi.

Materials and methods

Fungal strains and plant materials

A total of 11 strains belonging to four species were selected for pathogenicity analysis in May 2019 (Table 1). These strains were originally isolated from D. armandi in western China from July to August 2018 and May to July 2019 (Wang et al. 2022) and all D. armandi were collected from standing trees. All fungal strains used in this study were maintained in the China Forestry Culture Collection Center (CFCC, part of the National Infrastructure of Microbial Resources). Strains were grown on 2% malt extract agar (MEA; 20-g Biolab malt extract, 20-g Biolab agar and 1000-mL deionized water) in 6-cm plates for 10 days at 25 °C in the dark before being used for inoculation. The negative control was inoculated with 2% MEA without the fungi.

In total, 215 healthy 6-year-old P. armandii of a similar size (approximate height: 4–5 m and diameter: 7–8 cm) were selected from a nursery field in the Longcaoping Forestry Bureau of Foping County, Shaanxi Province (N 33°35'51", E 109°25'32"). In addition to considering the virulence of each species, a joint effect of different fungi that were frequently isolated from the same bark gallery or individual beetle was evaluated on the trees; therefore, both single-strain and two-strain concurrent inoculation methods (co-inoculation: L. qinlingense + O. shennongense) were used to evaluate the virulence of the ophiostomatoid strains. Leptographium qinlingense and O. shennongense are two dominant and widespread species associated with D. armandi; they are obtained from identical D. armandi or galleries of D. armandi. Therefore, we chose these two species for co-inoculation. Each of the 11 selected strains and a co-inoculation of the two strains from different species were inoculated into the stems of five trees.

Virulence tests

Each tree stem was subjected to the same treatment (inoculated with one strain or two co-inoculation strains) at the positions of three circling lines 35, 85 and 135 cm above the ground, and three inoculation holes were longitudinally distributed in each tree. Briefly, a 5-mm diameter hole was drilled horizontally into the phloem of P. armandii using a cork borer. Next, 5-mm mycelium plugs cut from 10-day-old actively growing colonies were inserted into the holes with a sterile toothpick (Yamaoka et al. 1998), allowing the contact of the mycelia with the xylem. The bark collected by the punch was used to cover the inoculation holes and was subsequently wrapped with parafilm and tape to prevent desiccation and air contamination.

Three batches of P. armandii saplings (among 215 P. armandii trees, 210 trees were divided into three batches, with 70 trees in each batch, for 11 single-strain inoculations, 2 co-inoculations and 2% MEA control treatments; another 5 healthy trees were not subjected to any inoculation and were used as healthy controls in the chemical analysis of host monoterpene compounds) that were subjected to the same treatment described above were used to evaluate the virulence of representative strains at 1, 30 (1 month) and 90 (3 months) days after inoculation (Fig. 1). At these time points, the phloems of the inoculated points were removed with a knife, and the length (mm) of necrosis on the cambium was measured and recorded.

Tissue samples from the peri-necrotic area were collected and placed in individual envelope bags for isolation of fungi under aseptic conditions. The tissue samples were cut into pieces of ~5 × 5 mm, disinfected using 1.5% sodium hypochlorite for 1 min, rinsed thrice with sterile water, and placed in 9-cm Petri dishes containing 2% MEA. All cultures were incubated in humid chambers in the dark at 25 °C and inspected daily for mycelial mass. The mycelium apex was used to purify all strains, and strains were identified by comparing culture characteristics.

Stain tests on xylem

Pathogens induce the host to produce oleoresin to block the ducts after fungal inoculation, which hampers water transport from the root to the crown. Thus, 3 months after fungal inoculation, the last batch of P. armandii saplings was cut at a circling line 5 cm above the ground and stained with a 5 mg/L of acid fuchsin solution for 24 h. The stems were then cut off at each inoculation line using an electric saw, the unstained part of the cross-section of the trunk was considered blocked (i.e. dry zone), and the blockage cross-sectional area was observed and photographed. Surface areas of the total cross-sectional and xylem blockages were determined using a computer-connected planimeter (ImageJ, Windows program contributed by George Silva).

Chemical analysis of host monoterpene compounds

Monoterpene volatiles were collected at 1 (31 May 2019), 30 (1 July 2019) and 90 days post-infection (d.p.i.) (1 September 2019) of P. armandii with fungi, a control and healthy trees; this was accomplished using a closed circulation dynamic headspace sampling system of Porapak-Q absorbent (i.e. 150 mg, 50–80 meshes in a Chrompack, 0.6 × 160 mm diameter long glass tube, CNW Technologies GmbH, Germany) through mini vacuum pumps (Atmospheric Sampling Instrument, QC-1S, Beijing Labor Protection Institute, China) at a 0.5-L/min air-flow rate for 30 min. The aeration sampling area was enclosed using a polyethylene (PE) film (23 × 20 cm, cut from a Reynold Oven bag, USA) and wrapped around the stem. Next, sampling tubes with volatiles trapped on the Porapak-Q absorbent were sealed at both ends with aluminum foil and placed inside a biological sampling box filled with dry ice (4 °C). The volatiles were extracted using 5 mL of HPLC-purified n-hexane. Piror to a gas chromatography (GC) analysis, the samples were concentrated to 40 μL using a weak nitrogen stream and stored at −20 °C. Three replicates were performed for each treatment group.

Samples were analyzed using gas chromatography-flame ionization detection (GC-FID; Agilent Technologies, Palo Alto, CA, USA) with an automatic sampler designed for liquid sample injections. Next, 1 μL of each aliquot was injected into an HP-5 column (Agilent Technologies, 30 m × 0.25 mm i.d. × 0.25-μm film thickness), with a nitrogen carrier gas flow at 1.8 mL/min at a 15-p.s.i. column head pressure. The flame ionization detector temperature was 270 °C, and the injector temperature was 250 °C. The oven temperature program started at 45 °C (isothermal, 3 min), increased linearly to 105 °C at 2 °C/min (isothermal, 1 min), and then to 250 °C at 15 °C/min (isothermal, 30 min). Absolute assays of monoterpenes were conducted using the following standards: S-(−)-α-pinene, camphene, β-pinene, myrcene, 3-carene, S-(−)-limonene and γ-terpinene, according to the retention time under the same conditions, and the monoterpene contents in the samples were calculated using an external standard.

Transcriptome analysis of P. armandii

Since all 6-year-old P. armandii were used for virulence tests of ophiostomatoid fungi, we selected 4-year-old P. armandii saplings to analyze the molecular changes in P. armandii after infection with a single inoculum or two concurrent inocula, including L. qinlingense (CFCC53937), O. shennongense (CFCC53921), L. qinlingense + O. shennongense (CFCC53937+ CFCC53921) and 2% MEA. All 4-year-old P. armandii stems underwent the same inoculation protocol at the position of the circling line 20 cm above the ground, and the inoculation method was used as described in the virulence tests (Fig. 1). A chemical analysis of the host monoterpene compounds revealed that the monoterpenoid release peaked at 30 d.p.i. Bark of P. armandii saplings inoculated with different inocula was subjected to transcriptome sequencing using the Illumina sequencing system at 30 d.p.i. Phloem tissue was collected from a 1–5 cm square at the disease/health junction, wrapped in aluminum foil, and placed in liquid nitrogen. These samples were stored at −80 °C in the laboratory prior to transcriptomic analysis.

Total RNA from the three biological replicate treatments was extracted using the RNAprep Pure Plant Kit (Polysaccharides and Polyphenolics-rich; TIANGEN, Beijing, China) according to the manufacturer’s protocol. RNA integrity and concentration were assessed using an Agilent 2100 Bioanalyzer (Agilent Technologies Inc., Santa Clara, CA, USA). The obtained RNA was stored at −80 °C.

The mRNA was isolated using the NEBNext Poly(A) mRNA Magnetic Isolation Module (NEB, E7490). The cDNA library was constructed using 1 μg from each sample following the manufacturer’s instructions for the NEBNext Ultra RNA Library Prep Kit (NEB, E7530) and NEBNext Multiplex Oligos (NEB, E7500) for Illumina. Briefly, the enriched mRNA was fragmented into ~200 nt RNA inserts, which were used to synthesize the first- and second-strand cDNA; double-stranded cDNA was subjected to end repair/dA-tail and adaptor ligation. Suitable fragments were isolated using Agencourt AMPure XP beads (Beckman Coulter, Inc.) and enriched by PCR amplification. Finally, the constructed cDNA libraries of the ophiostomatoid fungi were sequenced on a flow cell using an Illumina HiSeq™ sequencing platform. The entire set of raw reads was submitted to the Gene Expression Omnibus (GEO) at NCBI; the accession number is PRJNA1028765.

Low-quality reads (reads containing >5% unknown nucleotides or reads containing >50% of bases with a Q-value of ≤10%) were removed by developing a Perl script to obtain high-quality reads to ensure the accuracy of the subsequent transcriptome analyses. To generate nonredundant unigenes, clean reads were assembled de novo using the Trinity method with an optimized K-mer length of 25. Gene expression levels were estimated as fragments per kilobase of transcript per million fragments mapped using Cufflinks software.

DESeq and Q-values were employed to evaluate the differential gene expression (Love et al. 2014) using a model based on a negative binomial distribution. The false-discovery rate (FDR) control method was used to identify the P-value threshold in multiple tests to compute the significance of differences. In the present study, differentially expressed gene (DEG) parameters were set as follows: fold change ≥2 and FDR correction <0.01. For annotation, the unigenes were annotated based on the BLAST parameter: E-value ≤1e⁻⁵, HMMER parameter: E-value ≤1e⁻¹⁰. Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways were assigned to the assembled sequences using a Perl script. The KOBAS software (Mao et al. 2005) was used to test the statistical enrichment of DEGs in the KEGG pathways.

Statistical analysis

Column and scatter charts were generated using Origin 8.1, and data were analyzed using a one-way analysis of variance (ANOVA) (SPSS 19.0, IBM, Chicago, IL, USA). A least significant difference (LSD) test at the α = 0.05 level was used to evaluate the lesion lengths and blockage cross-sectional area of different strains/species during the same period. The lesion lengths of the same species at different time points (30 and 90 d.p.i.) were analyzed using Tukey’s test at the α = 0.05 level. Agilent ChemStation b.03.02 was used for data processing of monoterpenes to obtain the peak areas of major host volatiles. A standard curve established by the external standard transformed data for major host volatile contents was analyzed using a one-way ANOVA, followed by an LSD test at the α = 0.05 level. The relationship between lesion length, proportion of the blockage area, and total monoterpene content at 90 d.p.i. was analyzed using Pearson’s correlation. Transcriptomic analysis was performed using the BMK cloud platform, and the metabolic pathways used in this study were generated using Visio 2010.

Results

Pathogenicity tests

Healthy 6-year-old P. armandii were inoculated with 11 strains (Table 1) representing the following taxa: G. parakesiyea, L. qinlingense, O. shennongense and Ophiostoma sp. 1, according to 14 treatments, namely, a control (2% MEA), the single-strain inoculation of each strain and two-strain co-inoculations of O. shennongense + L. qinlingense (CFCC53932 + CFCC53941, CFCC53921 + CFCC53937). Three months after inoculation, no trees treated exhibited visible disease symptoms on the crown.

In the single-strain inoculation experiments, the necrotic lesions produced by L. qinlingense were significantly longer than those produced by the other species (Figs 2A, B and 3). The length of the lesion produced by Ophiostoma sp. 1 was shorter than that of L. qinlingense but longer than those of O. shennongense and G. parakesiyea. Furthermore, the length of the lesion produced by O. shennongense and G. parakesiyea did not differ significantly from that of the control (Fig. 2A and B). In the co-inoculation experiments of two strains/species, the lesion length induced by O. shennongense + L. qinlingense was significantly longer than that induced by the weakly pathogenic O. shennongense and shorter than that induced by the highly pathogenic L. qinlingense (Fig. 2A and B). There was no significant difference between the lesion lengths obtained after co-inoculation treatments with the two species (O. shennongense + L. qinlingense) (Fig. 2A).

Furthermore, there were no differences in necrotic lesion length between 30 and 90 d.p.i. for the control, G. parakesiyea, O. shennongense and Ophiostoma sp. 1 treatments; however, statistically significant differences were observed among the L. qinlingense and O. shennongense + L. qinlingense treatments (Fig. 2B).

For the weakly virulent G. parakesiyea and O. shennongense, the proportion of the blockage cross-sectional area was not significantly different from that of the control (Figs 2C and4). The highly virulent Ophiostoma sp. 1 and O. shennongense + L. qinlingense treatments had blockage cross-sectional area proportions that differed significantly from those of the control (Figs 2C and4). Additionally, the most virulent species among the six treatments, L. qinlingense, was associated with the largest proportion of blocked cross-sectional area, which differed significantly from the other five treatments (Figs 2C and4). So, the proportion of the blockage cross-sectional area waspositively correlated (R2 = 0.99, P < 0.01; Fig. 5) with fungalvirulence (lesion length).

Tissue samples were collected from the peri-necrotic area for isolation of the fungi under aseptic conditions. All four species of ophiostomatoid fungi that induced necrosis were successfully isolated from the necrotic area but not from the controls. Furthermore, O. shennongense and L. qinlingense, in which their co-inoculation induced necrosis, were successfully isolated from the necrotic area.

Analysis of host monoterpenes

As no significant differences in lesion length were observed following inoculation among the different strains of the same species, except for CFCC53941, the data generated from the different strains of the same species were chosen for two-strain treatments and were pooled for further analysis. A total of seven monoterpenes (S-[−]-α-pinene, camphene, β-pinene, myrcene, 3-carene, S-[−]-limonene and γ-terpinene) were analyzed from the n-hexane extracts. Overall, the total monoterpene content was positively correlated with fungal virulence (lesion length) at 90 d.p.i. (R² = 0.93, P < 0.01; Fig. 5). Compared with the healthy trees, the total contents of host monoterpenes, particularly O. shennongense, Ophiostoma sp. 1, L. qinlingens and O. shennongense + L. qinlingense treatments, initially increased from 1 to 30 d.p.i. and subsequently decreased up until 90 d.p.i.; the exception was the O. shennongense + L. qinlingense treatments, where the total monoterpene content continuously increased from 1 to 90 d.p.i. Notably, the total monoterpene content under each treatment, except for that in healthy trees and controls, was higher at 90 than at 1 d.p.i. (Fig. 6, Table S1 available as Supplementary data at Tree Physiology Online). In contrast to other monoterpenes, the γ-terpinene level peaked at 1 d.p.i. and decreased subsequently, and was detected only in healthy trees, the control and G. parakesiyea treatments at 3 months after treatment. Furthermore, 3-carene was not detected in the control, O. shennongense and Ophiostoma sp. 1 treatments at 30 and 90 d.p.i. (Table S1 available as Supplementary data at Tree Physiology Online). Among the seven monoterpenes detected from the artificial inoculation sites on the P. armandii trunks, the S-(−)-α-pinene levels were significantly high among control and healthy trees compared with the other six monoterpenes.

Regarding monoterpene contents at each measured time point (1, 30 and 90 d.p.i.) after inoculation, the L. qinlingense treatments produced the strongest response, and S-(−)-α-pinene, camphene, β-pinene, myrcene and S-(−)-limonene levels were significantly increased compared with those of healthy trees (Fig. 6, Table S1 available as Supplementary data at Tree Physiology Online). Furthermore, S-(−)-α-pinene levels were increased by >250 times compared with that of healthy trees at 30 d.p.i. (Table S1 available as Supplementary data at Tree Physiology Online). In addition, monoterpene content following Ophiostoma sp. 1 and O. shennongense + L. qinlingense treatments was significantly higher than that of healthy trees, the control, and G. parakesiyea treatment at 30 and 90 d.p.i. Similarly, the monoterpene contents of O. shennongense treatments were significantly higher than those of the healthy trees, control and G. parakesiyea treatments at 30 d.p.i. (Fig. 6). Collectively, these results indicated that among the five inoculation treatments, L. qinlingense and Ophiostoma sp. 1 induced a higher monoterpene accumulation than O. shennongense + L. qinlingense, O. shennongense or G. parakesiyea at 30 d.p.i. (Fig. 6, Table S1 available as Supplementary data at Tree Physiology Online). Although the content of all monoterpenes increased compared with the control and healthy treatments, other statistical data were not significant (Table S1 available as Supplementary data at Tree Physiology Online).

Global review of transcriptome sequencing data

Based on the finding that the co-inoculation-induced host produced a shorter lesion than the highly virulent fungi, 12 P. armandii samples (LQ: L. qinlingense treatments, OS: O. shennongense treatments, QS: O. shennongense + L. qinlingense treatments and control: 2% MEA treatments) collected following the L. qinlingense, O. shennongense, O. shennongense + L. qinlingense and 2% MEA inoculations were selected for RNA-seq; however, only nine P. armandii samples of LQ, OS and QS were used for comparative transcriptomic analysis (LQ vs QS and OS vs QS) (Table 2). A total of 57.69 Gb clean data with Q30 > 93.20% were obtained after adaptor removal and quality filtering; the clean data for each sample reached a 5.74 Gb size and above.

All clean reads of P. armandii were assembled into 54,556 unigenes with an N50 of 2157 bp; this indicated that the raw data demonstrated a high integrity for applicability in subsequent analyses (Table 2). The raw RNA-seq data were deposited in the Sequence Read Archive (SRA) of NCBI (SRA number: PRJNA1028765) under the BioSample accessions SAMN37852888–SAMN37852896.

All biological replicates in P. armandii samples were clustered using a hierarchical clustering analysis. Spearman’s Correlation Coefficient (r) was used as an evaluation index for correlation among biological replicates. The Pearson coefficient of all three biological replicates in all samples was >0.85, indicating that our sequencing data had high reliability and were suitable for further analysis (Fig. S1A available as Supplementary data at Tree Physiology Online). Hierarchical clustering analysis and principal component analysis confirmed that the replicates showed consistency, and checked the global patterns in the data, further verifying the quality of the dataset (Fig. S1 available as Supplementary data at Tree Physiology Online).

Differentially expressed genes in P. armandii

In the present study, 37,958 P. armandii functionally annotated genes were identified. The expression patterns of DEGs in the three groups (LQ, QS and OS) differed (Fig. 7A); OS was significantly similar to QS but different from LQ. To explore the differences in gene expression in P. armandii, the transcriptomes of the bark tissue in the LQ vs QS and OS vs QS groups were compared. In total, 269 DEGs were detected in the two comparisons, of which six were detected in the simultaneous comparison of LQ vs QS and OS vs QS (Fig. 7B). Compared with QS, 239 and 36 DEGs were detected in LQ- and OS-inoculated P. armandii, respectively, with 145 and 24 genes downregulated (Fig. 7C and D).

KEGG was used to predict the biochemical pathways of the DEGs; 62 signal pathways were categorized, and KEGG analysis enriched the DEGs in the top 10 metabolic pathways in the LQ vs QS group (P < 0.05) (Fig. 8 and Table S2 available as Supplementary data at Tree Physiology Online; Kanehisa et al. 2021). Among these, DEGs involved in flavonoid biosynthesis (ko00941), phenylpropanoid biosynthesis (ko00940), flavone and flavonol biosynthesis (ko00944), alpha-linolenic acid metabolism (ko00592), and stilbenoid, diarylheptanoid and gingerol biosynthesis (ko00945) were downregulated; in addition to ko00592, the other four biosynthetic pathways were related to the synthesis of phenolic compounds, which are influence host defenses against herbivores and pathogens.

Furthermore, nine signaling pathways were annotated, and only plant–pathogen interaction (ko04626) enriched the DEGs in the OS vs QS group (P < 0.05); DEGs were both upregulated and downregulated (Table S2 available as Supplementary data at Tree Physiology Online).

Differentially expressed genes involved in metabolic pathways

Based on the DEGs of KEGG, five DEGs in phenylpropanoid biosynthesis were downregulated in LQ vs. QS (−4.1 ≤ log2FC ≤ −1.463) (Fig. 9, Table S3 available as Supplementary data at Tree Physiology Online). Four enzymes, including phenylalanine ammonia-lyase (PAL), caffeic acid 3-O-methyltransferase/acetylserotonin O-methyltransferase (COMT), shikimate O-hydroxycinnamoyltransferase (HCT) and peroxidase (POD), showed downregulated expression. The downregulation of DEGs and enzymes primarily affected the biosynthesis of 5-hydroxyguaiacyl, p-hydroxyphenyl, guaiacyl and syringyl lignins.

A total of eight DEGs were associated with the flavonoid biosynthesis between the two groups (LQ vs QS); all the eight DEGs were downregulated (−11.093 ≤ log2FC ≤ −1.135) (Fig. 10, Table S3 available as Supplementary data at Tree Physiology Online). In both comparisons, enzymes critically associated with flavonoid biosynthesis included HCT, flavonoid 3′-monooxygenase (CYP75B1), bifunctional dihydroflavonol 4-reductase/flavanone 4-reduc-tase (DFR) and leucoanthocyanidin reductase (LAR). The downregulation of DEGs and enzymes primarily affected the biosynthesis of (+)-catechin and (+)-gallocatechin.

Furthermore, two genes associated with one enzyme (CYP75B1) were downregulated in flavone and flavonol biosynthesis, and two genes associated with one enzyme (HCT) were downregulated in diarylheptanoid and gingerol biosynthesis (Table S3 available as Supplementary data at Tree Physiology Online).

Three DEGs (one upregulated and two downregulated) were involved in plant–pathogen interactions in OS vs QS. The three genes belonged to three key enzymes: LRR receptor-like serine/threonine-protein kinase FLS2 (FLS2), LRR receptor-like serine/threonine-protein kinase EFR (EFR) and pathogenesis-related protein 1 (PR1) (Table S3 available as Supplementary data at Tree Physiology Online).

Discussion

Ophiostomatoid fungi associated with bark beetles are virulent to host trees and can exhaust tree defenses by inducing and utilizing the metabolism of terpenes, phenolics and other defense compounds (DiGuistini et al. 2011; Hammerbacher et al. 2013; Yamaoka 2017). We investigated the virulence of ophiostomatoid fungi associated with D. armandi infesting P. armandii and the resistance of the host to pathogenic fungi at the physiological and molecular levels.

Virulence trials of the four ophiostomatoid fungi were performed in healthy 6-year-old P. armandii. The results showed that all four species were virulent to P. armandii, and the level of virulence decreased in the following order: L. qinlingense > L. qinlingense + O. shennongense > Ophiostoma sp. 1 > O. shennongense > G. parakesiyea. Correlation analyses showed that the proportion of the blockage cross-sectional area and seven monoterpene contents in P. armandii were positively correlated with lesion length (Fig. 5). Furthermore, among the fungal communities associated with D. armandi-P. armandii, L. qinlingense exhibited the highest virulence (Figs 2–4). High virulence in this species has been reported previously. Following inoculation of L. qinlingense into P. armandii seedlings, the lesion length was 8.65 cm at 21 d.p.i. (Chen and Tang 2002; Pu et al. 2008; Hu et al. 2014; Thanh et al. 2014). Because L. qinlingense was detected in D. armandi infesting P. armandii at a relatively high frequency (Wang et al. 2022), it may have contributed to the damage observed in D. armandi-infected P. armandii.

A virulence analysis in the present study demonstrated that all four ophiostomatoid fungi blocked the xylem and inhibited water transport in the host (Fig. 4), with a higher virulence correlating with a higher proportion of blocked cross-sectional areas. Similar results were observed in 30-year-old P. armandii after inoculation with phytopathogenic ophiostomatoid fungi, which developed in the resin ducts; this effectively blocked the xylem resin ducts of P. armandii and concurrently produced blue-stained sapwood (Tang and Chen 1999; Chen and Tang 2002). Additionally, other pathogenic ophiostomatoid fungi can block the xylem of other coniferous trees (Wingfield et al. 2017). The dry zone of L. kaempferi (Lamb.) Carr reached the border between the sapwood and heartwood at 55 d.p.i. with Endoconidiophora fujiensis (Yamaoka et al. 1998), Endoconidiophora laricicola (Redfern et al. 1987) and G. clavigera (Solheim and Krokene 1998) generated larger necrotic lesions in the inner bark and dry zones of the sapwood of European larch and shore pines (Pinus contorta var. contorta); Endoconidiophora polonica penetrated much deeper into the sapwood than Ophiostoma penicillatum and O. aenigmaticum, and induced more extensive sapwood desiccation after inoculation Yezo spruce trees (Picea jezoensis) (Krokene and Solheim 1997; Yamaoka et al. 2000).

In the present study, the observed monoterpene levels from the time of fungal inoculation to 30 d.p.i. aligned with those of previous studies on conifers, in which the terpene content increased following inoculation with blue-stained fungi (Villari et al. 2012; Pan et al. 2018; Fang et al. 2020). Additionally, artificial inoculation trials have shown that total monoterpene and sesquiterpene concentrations in the phloem of the stem were significantly higher following inoculation of P. armandii with ophiostomatoid fungi (Ophiostoma sp., Leptographium sp. and L. qinlingense) (Tang and Chen 1999; Thanh et al. 2014), which is consistent with our experimental results. These results indicate that the fungi induced a host defense response. However, the energy of host trees is limited, and the trees face a trade-off in terms of allocation of limiting resources between growth and defense, which affects the ability to adjust to environmental challenges. Thus, there will be a marked decline in the induce defense in trees after an initial period, and secondary metabolite will altered along with it. The observed decrease in monoterpene content from 30 to 90 d.p.i. may reflect the changes of physiological responses to host defense induced by ophiostomatoid fungi on the temporal scales.

Additionally, the total monoterpene content of P. armandii changed based on fungal virulence. Highly virulent fungi (L. qinlingense and Ophiostoma sp. 1.) induced the host to released higher levels of monoterpenes than weakly virulent fungi (O. shennongense and G. parakesiyea) (Fig. 6 and Table S1 available as Supplementary data at Tree Physiology Online). These results are consistent with those of Pan et al. (2018), who found that the six predominant monoterpenes (α-pinene, camphene, β-pinene, myrcene, β-phellandrene and α-terpinolene) of Pinus yunnanensis, collected after exposure to highly virulent fungi, were significantly higher than those produced in response to weakly virulence fungi. Collectively, these results demonstrated that when different pathogens invade a host, the defense responses of a host depend on the nature of the invading pathogen; for example, highly virulent pathogens elicit stronger defense responses in comparison with pathogens with low virulence.

The results of co-inoculation showed that concurrent inoculation of P. armandii with L. qinlingense + O. shennongense reduced lesion length and total monoterpene levels compared with single L. qinlingense treatments but increased lesion length and total monoterpene levels compared with single O. shennongense treatments (Figs 2 and 6). In contrast, following a co-inoculation of Scots pine with nutritional (O. macrosporum) and pathogenic (O. clavatum) fungi, the lesion length under single-fungus inoculation was shorter than that of the co-inoculation treatment, and the local increase in (+)-α-pinene level was significantly lower than under co-inoculation (Villari 2012; Villari et al. 2012). There were synergetic effects of the two fungi on Scots pine. Besides, an inoculation with a mixture of E. polonica, O. penicillatum and Ophiostoma piceae showed almost the same results as an inoculation with E. polonica alone (Yamaoka et al. 2000). There are no synergetic or antagonistic effects of the three fungi on the trees; however, the mixture of L. qinlingense and O. shennongense showed antagonistic effects on P. armandii in the present study. Therefore, the transcriptome analysis (single inoculation vs co-inoculation) was performed to determine the regulation of gene expression related to disease resistance and immunity in P. armandii.

The response of P. armandii bark following fungal inoculation showed that the expression pattern of DEGs between QS and OS was similar, but varied from that of LQ (Fig. 7A). Most of the DEGs related to host defense in LQ vs QS were downregulated, especially DEGs associated with phenylpropanoid biosynthesis (ko00940) and flavonoid biosynthesis (ko00941); these DEGs are related to the synthesis of host defense substances, such as phenolic compounds (lignin and flavonoids). Moreover, although the accumulation of monoterpenes following inoculation with ophiostomatoid fungi increased (Fig. 6), the expression of their biosynthesis genes did not. This has been shown by metabolite and transcriptome analysis of Norway spruce bark, where monoterpenes and diterpene resin acids increased more rapidly after wounding in MeJA-treated than control; however, the expression of their biosynthesis genes did not (Mageroy et al. 2020a, 2020b).

The synthesis of secondary phenolic metabolites in plants is primarily accomplished through phenylpropanoid biosynthesis. Flavonoid biosynthesis is an important downstream pathway of phenylpropane biosynthesis (Li et al. 2019). In the comparison of LQ vs QS, seven enzymes (PAL, COMT, HCT, POD, CYP75B1, DFR and LAR), including 11 DEGs related to phenylpropanoid biosynthesis and flavonoid biosynthesis that encode the synthesis of p-hydroxyphenyl lignin, guaiacyl lignin, 5-hydroxyguaiacyl lignin, syringyl lignin, (+)-gallocatechin and (+)-catechin, were downregulated (Figs 9 and 10).

PAL is the primary and rate-limiting enzyme in phenylpropanoid biosynthesis and influences several metabolic pathways of secondary compounds, including flavonoids and lignin. PAL catalyzes and regulates the synthesis of several phenolic compounds and lignins; PAL1 and PAL2 are involved in the synthesis of lignin and flavonoids in Arabidopsis thaliana (Liang et al. 1980; Rohde et al. 2004). Pathogen stress can increase PAL activity; for instance, the expression of PAL is upregulated after infection with Rhizoctonia solani, which causes anaphylactic necrosis in maize leaves (Ma et al. 2015). Similar results were obtained in the present study; LQ treatments stimulated stronger expressions of PAL and DEG (BMK_Unigene_035696) than those of the QS treatments (Fig. 9).

Furthermore, several other key enzymes with downregulated expressions in LQ vs QS, also strongly influence plant resistance to pathogen invasion. For example, COMT and POD have integral functions in the pathways of lignin polymerization, pathogen resistance process, plant hormone metabolism and phenol catalytic oxidation process (Zubieta et al. 2001). HCT, a bifunctional enzyme, was simultaneously activated and expressed in phenylpropanoid biosynthesis and flavonoid biosynthesis (Figs 9 and 10).

Lignin acts as a natural physical barrier that limits pathogen invasion. In this study, the expression levels of key enzymes involved in lignin synthesis were higher in the LQ treatment group than in the QS treatment. The results of this study indicated that co-inoculation reduces the virulence of L. qinlingense from a molecular perspective. In addition, plants accumulate a large number of phenolic compounds to inhibit pathogens; this is associated with disease resistance (Buer et al. 2010; Panche et al. 2016). In this study, DEGs were downregulated in flavonoid biosynthesis in LQ vs QS, indicating that LQ treatments exerted greater pressure on plants than QS treatments, and more phenolic synthetic genes were activated to combat the invasion of L. qinlingense. This is consistent with the results of previous studies. After inoculation, E. polonica associated with I. typographus induced the expression of flavonoid substances in the host Picea abies and improved the induced resistance of the host to bark beetles and pathogens (Wadke et al. 2016). Ophiostoma bicolor inoculation of Picea koraiensis seedlings activated flavonoid biosynthesis in the host (Liu et al. 2022).

We did not detect lignin and phenolic substance levels; however, relevant studies have shown that the lignin level in the treated area was reduced in Scots pine after inoculation with the pathogen ophiostomatoid fungi (Villari 2012; Villari et al. 2012). Moreover, E. polonica can use phenolics as carbon sources for its growth to reduce the content of phenolic compounds (Hammerbacher et al. 2013; Wadke et al. 2016). In future studies, apart from molecular analysis, we will focus on the changes in these two substances in the host after pathogen inoculation.

In conclusion, in this study, the virulence of ophiostomatoid fungal species associated with D. armandi infesting P. armandii was verified at the physiological and molecular levels. Ophiostomatoid fungi can induce host defense responses and weaken host resistance, and the virulence of ophiostomatoid fungi may contribute to the invasion and colonization of P. armandii by D. armandi. Recently, the sustainable development of P. armandii in China has been increasingly threatened by D. armandi and associated ophiostomatoid fungi. This study provides valuable insights and a scientific theoretical basis managing of D. armandi and its associated fungi. Therefore, to determine whether ophiostomatoid fungi contribute to the invasion of conifers by bark beetles and evaluate the consistency of the results obtained, additional bark beetle–fungal pathogenicity studies are warranted.

Acknowledgments

We thank Professor Zhongdong Yu from Northwest A&F University, Mr Jiaxi Yi from Sheng Nong Jia forestry pest natural enemy breeding farm, Hubei, and Mr Anmin Li from Xiaolong Mountain Forest Experiment Bureau of Gansu Province, and Tianshui, for their assistance with insect sample collection.

Authors' contributions

H.W. designed the study, analyzed the data, contributed to experiment design and wrote the manuscript. H.W. and T.W. collected the samples. H.W., Y.L. and D.L. performed DNA extraction and PCR amplification. Q.L. reviewed and approved the final manuscript. All authors have read and agreed to the published version of the manuscript. All authors contributed to the article and approved the submitted version.

Funding

National Natural Science Foundation of China (project no. 32230071, 32071769, 32201569).

Conflict of interest

The authors declare no conflict of interest.

Data availability

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary data.

References