1. Introduction
Scots pine (
Pinus sylvestris L.) is the most widely distributed pine species and is found throughout Eurasia [
1]. It is a species of major economic importance, widely used in timber, pulp, and paper production [
2]. It is also a keystone species providing stability for large ecosystems [
3,
4] and has been a dominant species for millennia [
5].
One of the main causes of mortality of pine, excluding bark beetles – genera Ips De Geer, 1775 and Tomicus Latreille, 1802 [
6,
7], are Heterobasidion complex fungi, in particular, the basidiomycete Heterobasidion annosum (Fr.) Bref., which causes root rot [
8]. In the European forest sector, the Heterobasidion complex causes combined economic losses in the hundreds of millions of euros annually [
9]. While the exact proportion of the losses specifically related to Scots pine is not known, given the widespread distribution of Scots pine in European forests and the high economic value of Scots pine, the economic impact of this pathosystem on Scots pine is significant.
Currently, mitigation of
Heterobasidion root rot includes forest management activities (monitoring, creation of sustainable forests, sanitary felling, etc.) [
10,
11,
12,
13]. One of the strategies for the protection of renewed pine stands is the use of biological control agents based on the fungus
Phlebiopsis gigantea (Fr.) Jülich [
14,
15,
16]. Breeding for resistance or tolerance is promising, as genetic components for
H. annosum resistance have been determined [
17], Scots pine clones with varying resistance against
H. annosum have been described [
18] and the heritability and genetic gain values for breeding for resistance against root rot have been calculated [
19].
To increase the efficiency of breeding programs, genetic mechanisms and candidate genes for resistance need to be identified, thus, studies of gene expression changes in the host in response to inoculation and related genetic testing have been performed for Scots pine reacting to
H. annosum, and also for related pathosystems [
20,
21,
22,
23,
24]. Genetic virulence factors of
H. annosum have also been studied [
25,
26]. A study combining genomic and transcriptome approaches identified genes linked with secreted proteins as potential virulence factors (besides genes involved in oxidation- reduction processes and genes encoding domains relevant to transcription factors) [
27]. In the same study secretome annotation and analysis in the pathogen-host interactions database [
28] showed that the most virulent
H. parviporum isolate was found to contain many carbohydrate-active enzyme genes for cell wall degradation and an increased amount of secreted proteins for lignin degradation. Proteins of the reactive oxygen species-scavenging system and proteases were identified as an important part of the secretome. This study also identified cytochrome P450 proteins as significant for pathogenesis, which could be explained by the involvement of cytochrome P450 gene family members in detoxification of substances from the environment and synthesis of fungal toxins [
29]. However, the study of Zeng et al. (2018) did not identify a specific gene as the sole determinant for variations in virulence between different
H. parviporum isolates.
Recently, dual transcriptome studies have provided information on host and pathogen transcriptomes simultaneously [
30,
31,
32,
33]. However, not all of these studies focus on the host-pathogen interaction [
31]. The dual transcriptome studies are also possible only if high number of reads can be obtained from samples to provide sufficient amount of pathogen transcriptome reads, as the proportion of transcripts from the pathogen can be below 1% [
32]. No dual transcriptome studies on the
H. annosum –
P. sylvestris pathosystem have been published and, to best of our knowledge, there are no investigations providing data on the
H. annosum (
sensu stricto) transcriptome during infection of Scots pine. However, in a study by Lunden et al., 2015 about inoculation of Norway spruce with
H. annosum s.s., delta-12 fatty acid desaturase and clavaminate synthase were indicated as potential virulence factors [
33]. Polysaccharide-degrading enzymes and lignin-degrading enzymes have long been regarded as important for pathogenesis [
34]. Expression of these genes was also detected in the dual transcriptome experiment by Kovalchuk et al., 2019 [
30]. The same study also mentions that seven
Heterobasidion genes classified as lipases are expressed, but no further details were provided in that study. Expression of lipid metabolism genes (encoding lipases, malic enzyme) was detected in
Fusarium circinatum Nirenberg & O'Donnell while infecting
Pinus species [
32]. In addition to expression of polysaccharide- and lignin-degrading enzyme genes, accumulation of toxins and oxalate have been reported in host tissues colonized by
Heterobasidion [
34].
Previous research has shown that analysis of fungal transcriptomes enables understanding of adaptations by pathogens to different host species [
32]. Therefore, accumulation of
in planta gene expression data for
H. annosum can enable analysis of the variation of defense strategies within-the host species by inoculation experiments with genetically characterized fungal isolates. This could be used as a basis for the identification of molecular markers for resistance-oriented forest tree breeding. In addition, comparison of fungal gene expression between isolates could provide insights into differences in aggressiveness or virulence, providing a better understanding of the role of genetic variation of pathogens in the development of plant diseases. This investigation provides an initial assessment of time post inoculation-dependent differentially expressed pathogen genes from one to four weeks post inoculation. This assessment is based on mapping a large amount of reads to the
H. annosum transcriptome (> 11 million). This provides information about genes and pathways important in early pathogenesis, and demonstrates the feasibility of the utilized sequencing technology.
2. Materials and Methods
One year old Scots pine saplings were inoculated with a Latvian isolate (ID HA2) of
H. annosum s.s.. The saplings were obtained from the Latvian State Forests seedling nursery in Kalsnava, Latvia, and were grown from mixed origin improved seed material. For inoculation, a wooden plug grown-through with the
H. annosum isolate was placed on an area of the stem with the bark surface removed. Samples for RNA extraction were collected at one, two, three, and four weeks post inoculation and stored at -80 °C until RNA extraction. Four biological replicate samples were collected for each time point. RNA extraction was performed following a CTAB-based protocol [
35]. The sampling site included the inoculation site. These time points were chosen, firstly, to allow time for the pathogen to grow into the host tissue and form biomass to increase the amount of obtainable reads and, secondly, to allow partial comparison with previous studies (Lundén et al., 2015 [
33] and Zamora-Ballesteros et al., 2021 [
32]) which used 5 and 4 days post inoculation as their only time points, respectively. In context of infection progression, in the related
Picea abies (L.) H. Karst. –
H. parviporum Niemelä & Korhonen pathosystem, germ tubes form in roots within 24 h after spore adhesion, colonization of cortical tissue happens 24 – 48 post inoculation, and the endodermis is reached at 72 h [
36]. By 12 to 15 days post inoculation stelar cells are deteriorated but plant responses involving papillae formation and lignification in cortical and endodermis tissue in roots can be observed [
36]. Thus, the selected time points represent early infection and later stages.
RNA extraction was followed by quality control of the RNA samples using the Agilent 2100 Bioanalyzer and RNA 6000 nano kit (Agilent, Cat. No. 5067-1511) for RNA integrity determination and quantitation using the Qubit spectrofluorometer and Qubit RNA BR Assay Kit (Thermo Fisher Scientific, Cat. No. Q10210). The obtained material was quite degraded. RIN values ranged from 1.9 to 4.7 (average 2.6) and three samples didn’t produce a RIN value but were included as good results were obtained using qPCR to detect the presence of
H. annosum RNA in the total RNA samples. Considering that short read (100bp) sequencing was used, thus reducing the influence of degraded RNA, sequencing libraries were prepared with all samples. RNA quality control data are provided in
Supplementary Table S1.
To confirm presence of
H. annosum RNA in the extracted RNA samples, one-step RT-qPCR with
H. annosum-specific primers targeting the laccase gene was performed. The GoTaq® 1-Step RT-qPCR System (Promega, Cat. No. A6021) was used with the forward primer: 5'-CCAGAAAGTAGACAATTATTGGATTCG-3', and reverse primer: 5’-GAGTTGCGGCCATTATCGA-3’ [
37]. Reaction mixture per sample was: 10 µl GoTaq qPCR Master Mix (2X); 0.4 µl GoScript RT Mix for 1-Step RT-qPCR (50X); 0.33 µl CXR Reference Dye; final concentration of each primer 200 nM; nuclease-free water to a volume of 20 µl. Thermal cycling profile: 37 °C 15 min, 95 °C 10 min, followed by 40 cycles of 95 °C 10 s, 60 °C 30 s, 72 °C 30 s, then followed by a melt curve stage during which the temperature is increased from 60 °C to 95 °C in 0.3 °C intervals. One-step RT-qPCR was performed in an Applied Biosystems StepOnePlus qPCR machine. To avoid the influence of primer dimers on the results, melting curve peak height for the specific qPCR product was used instead of Ct value. This approach was possible because none of the reactions reached plateau phase. Peak height of the specific product was divided by RNA concentration to obtain a ratio characterizing the proportion of
H. annosum RNA in the RNA sample.
Transcriptome sequencing libraries were created using the MGIEasy RNA Directional Library Prep Set (MGI Tech Co., Cat. No. 1000006385) following the manufacturer’s protocol. 500 ng of total RNA was used for each library. Ribosomal RNA was depleted by use of the Qiagen QIAseq FastSelect –rRNA Plant Kit (Cat. No. 334315). Incubation was performed as described in
Table 4 of the QIAseq® FastSelect™ Handbook (according to treatment of RIN < 3 samples). The reaction mixture for this step was as following: 10 µl of RNA sample containing 500 ng of RNA, 4 µl fragmentation buffer (MGI), 4 µl directional RT Buffer 1 (MGI), 1 µl diluted directional RT buffer 2 (MGI), 1 µl QIAseq FastSelect −rRNA Plant reagent. RT enzyme mix (MGI) was only added after this incubation, before the reverse transcription step. The rest of the library preparation protocol was not modified from the manufacturers’ protocol. DNBSEQ-G400RS High-throughput Sequencing Set (MGI Tech Co., Cat. No. 1000016950) was used for sequencing. Sequencing was performed on a DNBSEQ G400 sequencer (MGI Tech Co.) by the Latvian Biomedical Research and Study Centre (Riga, Latvia) on May 20, 2023. The manufacturer’s standard protocol was used for sequencing.
The CLC Genomics Workbench v.23 was used for data import and CLC Genomics Workbench v.21.0.5 (including the Blast2GO commercial plugin v.1.21.14 (BioBam Bioinformatics S.L.)) was used for data analysis and annotation. Data analysis workflow:
demultiplexing (automatically performed by the sequencer),
trimming (settings: quality limit: 0.05; max. number of ambiguous nucleotides: 2; no adapter trimming; trim homopolymers; no terminal nucleotide removal; max. Length: 150),
RNA-seq analysis (settings: no spike-in controls; mismatch cost: 2; insertion cost: 3; deletion cost: 3; length fraction: 0.8; similarity fraction: 0.8; no global alignment; strand specific: both; library type: bulk; max. number of hits for a read: 10; count paired reads as two: no; expression value: total counts; create reads track: yes),
differential gene expression analysis (settings: technology – whole transcriptome RNA-Seq; filter on average expression for FDR correction: No; metadata table: yes; test differential expression due to: WPI (weeks post inoculation); while controlling for: not set; comparisons: all group pairs), selection of the differentially expressed genes to perform
annotation using the Blast2GO plugin (settings: blast program: blastx-fast; Blast DB = nr; mapping using Goa version 2022_08). Mapping of reads to a reference transcriptome of
H. annosum was performed [
38,
39] (settings: masking mode – no masking; update contigs: no; match score: 1; mismatch cost: 2; cost of insertions and deletions: linear gap cost; insertion cost: 3; deletion cost: 3; length fraction 0.5; similarity fraction: 0.8; no global alignment; auto-detect paired distances; non-specific match handling: map randomly) to be able to calculate the percentage of
H. annosum reads. The reads were also mapped against the transcriptome of
P. sylvestris from Wachowiak et al., 2015 [
40] using the same settings.
A gene (transcript) was considered differentially expressed if the comparison between time points produced a statistically significantly different count of reads mapped to the transcript in question with the p value below 0.01. Results are expressed as fold change comparing different time points (with four biological replicates for each time point). Calculations of the differential gene expression related parameters were performed using the differential gene expression analysis tool in the CLC genomic workbench software, see paragraphs 31.6 and 31.6.4 in the manual [
41] for more details.
For visualization of the differential expression, the online tool for up to 6-group Venn diagram creation InteractiVenn [
42]. Visualization of other data was performed using the tools in CLC software and Blast2GO plugin.
4. Discussion
Differences among the libraries in the proportion of H. annosum reads might indicate biologically relevant information such as the proportion of living pathogen in the tissue from which RNA was extracted but other explanations like sampling effects, host genotype effects, unequal rRNA depletion or induced host rRNA production at some time points could influence the results. This could also be an indication of a varying positive effect from the inoculum plug still providing resources to the pathogen but this hypothesis was not explicitly analyzed in this study. If one outlier (library with largest difference from average percentage of H. annosum reads per group) is removed per time point, a strong correlation (R2 = 0.8116) between longer time since inoculation and lower proportion of H. annosum reads is observed. Yet, given the many factors influencing the proportion of H. annosum reads, we refrain from conclusions.
Most of the differentially expressed genes lack detailed annotations, however, some of the genes showing the highest upregulation during infection might have a role in disturbing host cell integrity (e.g. carotenoid ester lipase precursor, which is secreted and can cross the cell wall, [
43], and erylysin B [
44]), or can affect lipid metabolism (e.g.carotenoid ester lipase precursor, fatty acid desaturase-domain-containing protein, delta-12 fatty acid desaturase [
33], elongase of fatty acids ELO, malic enzyme [
45]), terpene metabolism (terpenoid cyclases/protein prenyltransferase alpha-alpha toroid), are involved in oxidation-reduction processes (aldo/keto reductase [
46]), cellular growth (GPI mannosyltransferase 3 [
47]), cell development and metabolism control (methionine adenosyltransferase (synonym for S-adenosylmethionine synthetase) [
48]) and stress tolerance (heat shock protein 70). Increased expression of delta-12 fatty acid desaturases and related protein genes is especially pronounced comparing 1 WPI to 4 WPI. This could be related to lignin degradation [
49], especially in the context of elevated transcription of polysaccharide lyase [
50]. The initial stage of infection requires the fungal hyphae to penetrate host cell walls to colonize the tissue. As the structural integrity of conifer tissue is mainly ensured by cellulose, hemicellulose and lignin, polysaccharide-degrading enzymes are needed to penetrate the tissue and to access nutrition sources away from potential competitors on tissue surface. Several microorganisms in soil and tissue surfaces can negatively affect the growth of
H. annosum, therefore the ability to quickly penetrate cell walls and to colonize tissues is important for survival of the fungus and successful infection of the host [
34].
At one week post inoculation, a terpenoid cyclase gene is more highly expressed than at two or three weeks post inoculation. However, as terpenoids represent the most diverse and abundant class of natural products and have high functional diversity[
51,
52], any suggestions about the possible role in the infection process would be highly speculative.
Genes downregulated during early stage infection have roles in transcription and translation (transcription regulator, 40S ribosomal protein S26, RS27A protein), polyamine regulation (ornithine decarboxylase antizyme-domain-containing protein [
53]), protein folding as chaperones (groes-like protein, HSP20-like chaperone), cell development / signaling (Cell division control/GTP binding protein). Downregulation of the negative regulator of differentiation 1 one week after inoculation could result in positive regulation of differentiation.
Clearly, improved characterization of the proteins encoded by the differentially expressed genes is needed to obtain a better understanding of the molecular processes affecting
H. annosum pathogenicity. However, the results obtained about annotated genes suggest that delta-12 fatty acid desaturases could be a virulence factor. The delta-12 fatty acid desaturase gene expression in
H. annosum – Norway spruce pathosystem were detected by Lundén et al. (2015) [
33]. These authors suggest a central role of this enzyme in sustaining of fungal growth. Considering the extremely small amount of
H. annosum transcripts in the work of Lundén et al. and the small proportion of
H. annosum reads in our work, identification of delta-12 fatty acid desaturase gene as important for the pathogen indicates that this gene plays an important role during infection. An increase in the expression of other transcripts influencing lipid metabolism was also observed (malic enzyme, elongase of fatty acids and an unspecified type of fatty acid desaturase-domain-containing protein). The gene for malic enzyme was the only differentially up-regulated gene comparing the one week post inoculation time point against all other time points. This concurs with findings showing the role of fungal lipids and lipid biosynthesis proteins in plant-pathogen interactions as potential virulence factors [
54,
55]. Malic enzyme is potentially highly interesting in this context. As reviewed by Voreapreeda et al., (2013) [
56], there are different types of malic enzymes (EC 1.1.1.38 – EC 1.1.1.40), but the malic enzyme gene identified in our study (sequence ID: CCPC2187.b1 [
39]) was not categorized to a specific group. However, activity of malic enzyme can drastically change the lipid content of oleaginous fungi [
57]. Furthermore, a study investigating the
H. annosum transcriptome during saprotrophic growth on Scots pine bark, sapwood and heartwood detected down-regulation of malic enzyme during growth on bark and heartwood chips [
58]. This may suggest that malic enzyme activity is needed specifically when invading a living host, however further research is needed, as in the previous study, the time post treatment was 3 months, in contrast to one week in our study.
The observed downregulation of glycoside- and glycosyl hydrolase genes at 1 WPI compared to 4 WPI might seem contra-intuitive, as enzymes of these groups and carbohydrate-active enzymes in general have been recognized as one of the main components facilitating fungal invasion [
59,
60]. Glycoside- and glycosyl hydrolases have a role in degradation of plant cell wall components, thus it could be expected to see increased activity of these genes at earlier stages of infection. However, this class of proteins is vast and complex [
61], and therefore these could be family members with other functions. Members of these groups could be involved in fungal cell wall remodeling, putting it in context with increased fungal growth at later stages of infection.
Up-regulation of specific transcription factors sensitive to fungal or plant hormones was not observed, except for one transcription regulator at 1 WPI with limited information about this protein. In general, the biological functions of most highly expressed genes specific to 1 WPI indicates attempts to mitigate oxidative stress and other elements involved in plant defenses. Response to oxidative stress could be linked with protection from reactive oxygen species involved in host defense responses [
62]. Alternatively, it could also be associated in the pathogen’s own development-related cell signaling [
63] or both. Other biological processes of 1WPI-specific transcripts are related to metabolism, gene expression and response to heat, osmotic and oxidative stress. In addition, the glyoxylate cycle can also facilitate reduction of oxidative stress [
64]. In a broad sense, these transcriptome changes suggest that the pathogen is employing protective measures against plant defenses.
Considering the high proportion of differentially expressed genes lacking detailed annotations, there are still many unknown factors and genes that are involved in the pathogenicity mechanism. More research on the molecular biology of H. annosum is needed to gain better understanding of virulence factors. Further research using genetically identical host plants to obtain a high-quality dual transcriptome of H. annosum – P. sylvestris interaction at different time points post inoculation can provide additional information about the interactions between the pathogen and host during the infection process.
Author Contributions
Conceptualization, M.R., D.R. and V.Š.; methodology, M.R., D.R. and V.Š.; software, M.R., and V.Š.; validation, M.R., and V.Š.; formal analysis, M.R., and V.Š.; investigation, M.R., and V.Š.; resources, M.R., and V.Š.; data curation, M.R., and V.Š.; writing—original draft preparation, M.R., and V.Š.; writing—review and editing, M.R., D.R. and V.Š.; visualization, M.R., and V.Š.; supervision, D.R. and V.Š; project administration, M.R., and V.Š.; funding acquisition, M.R., and V.Š. All authors have read and agreed to the published version of the manuscript.