Skip to main content
  • Research article
  • Open access
  • Published:

The transcriptome of Pinus pinaster under Fusarium circinatum challenge



Fusarium circinatum, the causal agent of pitch canker disease, poses a serious threat to several Pinus species affecting plantations and nurseries. Although Pinus pinaster has shown moderate resistance to F. circinatum, the molecular mechanisms of defense in this host are still unknown. Phytohormones produced by the plant and by the pathogen are known to play a crucial role in determining the outcome of plant-pathogen interactions. Therefore, the aim of this study was to determine the role of phytohormones in F. circinatum virulence, that compromise host resistance.


A high quality P. pinaster de novo transcriptome assembly was generated, represented by 24,375 sequences from which 17,593 were full length genes, and utilized to determine the expression profiles of both organisms during the infection process at 3, 5 and 10 days post-inoculation using a dual RNA-sequencing approach. The moderate resistance shown by Pinus pinaster at the early time points may be explained by the expression profiles pertaining to early recognition of the pathogen, the induction of pathogenesis-related proteins and the activation of complex phytohormone signaling pathways that involves crosstalk between salicylic acid, jasmonic acid, ethylene and possibly auxins. Moreover, the expression of F. circinatum genes related to hormone biosynthesis suggests manipulation of the host phytohormone balance to its own benefit.


We hypothesize three key steps of host manipulation: perturbing ethylene homeostasis by fungal expression of genes related to ethylene biosynthesis, blocking jasmonic acid signaling by coronatine insensitive 1 (COI1) suppression, and preventing salicylic acid biosynthesis from the chorismate pathway by the synthesis of isochorismatase family hydrolase (ICSH) genes. These results warrant further testing in F. circinatum mutants to confirm the mechanism behind perturbing host phytohormone homeostasis.


Plant hormones play a crucial role in plant-pathogen interactions, especially salicylic acid (SA), jasmonic acid (JA) and ethylene (ET) which are primary signals for induction of defense response. Generally, ET and JA act synergistically in response to the infection of necrotrophic fungi, while SA is most commonly expressed in response to biotrophic or hemibiotrophic fungi [41]. However, this is only a simplistic classification and crosstalk between phytohormones are much more complex and context dependent. Although an antagonistic relation between JA and SA has been reported [59, 76, 95], this antagonism seems to be absent in the defense response of Arabidopsis to Plectosphaerella cucumerina [8], Pseudomonas syringae and Peronospora parasitica [23]. A cooperation between the two phytohormones has also been described in Picea abies [4] and Zea mays [36]. A role of auxins in plant-pathogen interaction has also been reported, modulating signaling pathways of other hormones, resulting in positive or negative effect in resistance [5, 58, 73, 77]. Considerable effort has been dedicated to understanding phytohormone signaling in plant defense, while knowledge on the role of fungal hormone production is limited.

A remarkable aspect of Fusarium species in the Fusarium fujikuroi species complex (FFC), is the ability to synthesize phytohormones, including gibberellins [12, 101] and auxins [103], that contribute to plant disease. However, the underlying molecular mechanism as well as their role in plant interactions remains unclear. Two mechanisms have been suggested: perturbing plant processes to favor invasion and nutrient uptake, and/or acting as signals for the fungus to engage appropriate physiological processes to allow adaption to the new environment [21]. Gibberellic acid (GA) production has been well described in the rice-infecting fungus Fusarium fujikuroi [12] and a correlation between GA levels and virulence has been reported [29]. GA biosynthetic genes are organized in a gene cluster, and while most of the species of the F. fujikuroi species complex have the entire GA biosynthetic gene cluster, F. circinatum was reported to have only one gene [12, 64]. Indol-3-acetic acid (IAA), the most common form of auxins, can be synthesized from tryptophan by the indol-3-acetamide (IAM) pathway, and IAM-related genes have been reported in four Fusarium pathogenic fungi: Fusarium verticillioides, Fusarium oxysporum, F. fujikuroi and Fusarium proliferatum [103]. In the same study, the deletion of an IAM-related gene resulted in drastic reduction of IAA production in F. proliferatum. Similarly, F. oxysporum transgenic lines containing two IAM genes produced significantly more IAA than the wild type when infecting Orobanche, leading to enhanced virulence [24]. ET producing fungi range from necrotrophic, like Botrytis cinerea, to biotrophic, such as F. oxysoporum f. sp. pini, and a role in perturbing the plant phytohormone homeostasis has been suggested [21]. Silicon treatment has been shown to induce brown spot resistance in rice against Cochliobolus miyabeanus by disarming fungal ET [105]. In the case of Colletotrichum sp. pathogens, ET is required for the formation of appressoria [38]. To our knowledge, expression of F. circinatum genes related to hormone biosynthesis or signaling, besides GA, has not been studied.

Fusarium circinatum Nirenberg & O’Donnell is described as one of the most important pathogens worldwide, affecting more than 60 Pinus species as well as Pseudotsuga menziesii Mirb. (Franco) [118]. The fungus can cause damage to seedlings as a wilt disease but also to mature trees, known as pitch canker disease. In Europe, the pathogen is currently present in the Atlantic area of Spain [87] and in Portugal [14], where Pinus species are grown. The maritime pine (Pinus pinaster Ait.) is the dominant species in the Mediterranean area with more than 2.3 million hectares [84], especially in the Atlantic coast of France, Portugal and Spain. Although P. pinaster has shown moderate resistance to the pathogen (mean lesion length of 5 mm compared to 28 mm in P. radiata [49] and a 66% incidence rate in a provenance/progeny trial of artificially inoculated seedlings [35], the presence of the fungus in the area represents a serious threat to this species. In spite of the special effort done in management of pitch canker disease based on cultural, biochemical and biological control, the pathogen has still not been eradicated. Therefore, selection of genetically resistant genotypes against F. circinatum seems to be an appropriate approach for disease management.

Advances in sequencing approaches have allowed the generation of new resources for conifers. RNA sequencing (RNAseq) approach allows the characterization of the transcriptome even in species for which no reference genome is available or is incomplete. In both instances, RNAseq reads can be assembled de novo into a transcriptome [17, 83, 110]. A variant of this technology is the dual RNAseq, which captures the transcriptome of the host and pathogen simultaneously, so that gene expression in both organisms may be determined [45, 66, 72, 116].

Various transcriptional changes are apparent in the Fusarium circinatum-Pinus pathosystem. Transcriptome analysis of Pinus radiata inoculated seedlings was recently published [18] showing induction of genes related to abscisic acid (ABA) signaling, auxin responsive-like proteins, gibberellin-regulated protein precursor, as well as induction of pathogenesis-related (PR) proteins, phosphorylase family protein (PFP) and genes related to physical and chemical barriers to restrict pathogen invasion. Similarly, Donoso et al., [32] detected an up-regulation of genes encoding thaumatin-like protein (PR5) in P. radiata inoculated seedlings by a RT-qPCR assay. Davis et al., [27] reported SA and JA induction of chitinase (PR3) in pine seedlings inoculated with Fusarium subglutinans f. sp. pini, suggesting a potential role for PR proteins in pine defense.

While most pine studies have focused their attention on host defense [18, 32], little is known about F. circinatum genes involved in pathogenicity. Muñoz-Adalia et al. [70] suggested 5 candidate genes that could be involved in F. circinatum virulence based on high similarity with other Fusarium species. Recently, using a dual RNAseq approach, Visser et al., [109] reported differences in the expression of F. circinatum ergosterol biosynthesis genes during the infection of Pinus tecunumanii and Pinus patula seedlings, suggesting a role for this pathway in pathogen virulence. The study also alluded to a role of phytohormone signaling in pine defense.

The aim of this study was to elucidate the role of phytohormones in moderate resistance of P. pinaster to F. circinatum and determine key steps where the pathogen could be manipulating host phytohormone balance to its own benefit, leading to host susceptibility. For this purpose, we determined by a dual RNAseq approach the expression profile of both organisms during the host-pathogen interaction at different times after inoculation (3, 5 and 10 days post-inoculation, dpi). Furthermore, a good quality P. pinaster de novo transcriptome assembly was generated, improving current P. pinaster genetic resources.


Pathogen colonization and symptom development

Fusarium circinatum was observed growing in the resin drop of P. pinaster seedlings during the first three days. At 4 dpi the pathogen had entered host tissue in only one out of six plants analyzed while at 5 dpi, the pathogen had penetrated all of them. During the following days, the fungus continued growing within the host and progressed from the tip along the stem. At 8 dpi some of the plants showed visible damage at the inoculation site. At 11 dpi inoculated seedlings had a lesion length of approximately 5 mm from the tip. The sampling times chosen for the RNAseq assay were: 3 dpi (the fungus had not penetrated within the host and was growing in the resin drop), 5 dpi (the fungus penetrated within host tissue) and 10 dpi (lesion was visible in the shoot tip of all seedlings). Figure 1 shows example images from these time points and the progression of the fungus over the different time points.

Fig. 1
figure 1

Fusarium circinatum inoculation of Pinus pinaster. a-c: symptoms at the shoot tip of inoculated Pinus pinaster seedlings at each sampling time point. a: 3 dpi, no visible symptoms. b: 5 dpi, no visible symptoms. c: 10 dpi, symptoms at the inoculation site. e, g, i: transverse sections of P. pinaster inoculated shoot tip visualized under epifluorescence microscope at each sampling time after inoculation. d, f, h: growth of the fungus in the resin drop at each time point. d-e: 3 dpi. f-g: 5 dpi. h-i: 10 dpi. White arrows indicate colonization of the fungus within host tissue. Bar 100 μm. j: progression of a green fluorescent protein (GFP)-tagged Fusarium circinatum strain within Pinus pinaster shoot tissue at different days after inoculation

By the end of the experiment, all inoculated seedlings showed symptoms of disease including discolored brown stems, necrosis and needle desiccation at the tip (Additional file 1). Inoculated seedlings showed a lesion length at the tip of 1.5 cm ± 0.59 (SD) and F. circinatum was re-isolated from all 6 tips cultured on Fusarium Specific Media. Mock-inoculated seedlings did not show symptoms of disease and the fungus was not recovered from any of them.

Pinus pinaster de novo reference transcriptome

Quality of preliminary assemblies

A total of 21 different Trinity and TransABySS preliminary assemblies were built, showing differences in quality based on the parameters used (Additional file 2). None of the assemblies produced transcripts with unknown bases and the GC content were similar between them. In silico normalization of Trinity assemblies produced the highest number of transcripts, mean length and better N50 values. However, in non-normalized assemblies more fragments were mapped and the percentage of good mapped contigs was higher (Additional file 3). Therefore, all Trinity assemblies (normalized and non-normalized) were used for building the final de novo transcriptome assembly. When comparing Trinity and TransABySS assemblies at the same kmer value, Trinity showed better quality statistics, therefore TransABySS assemblies with the same kmer value as Trinity were discarded. As Trinity does not allow the use of kmer values higher than 32, all TransABySS assemblies with higher kmer values were conserved.

A total of 18 good quality preliminary assemblies were used as input for the EviGene pipeline to build the P. pinaster de novo transcriptome. When merging best quality assemblies, EviGene pipeline classified 90.7% of the sequences as redundant and uninformative, which means that 49,624 sequences (9.3%) were non redundant coding genes and were used to generate the P. pinaster transcriptome (Additional file 2).

Benchmarking Universal Single Copy Ortholog (BUSCO) analysis against the embryophyta_odb9 lineage database identified 1261 (87.57%) complete BUSCOs in a total of 1440 BUSCO groups searched, from which 1109 (77.01%) were single copy, 152 (10.56%) duplicated, 145 (10.07%) missing and 34 (2.36%) fragmented. Regarding the eukaryote_odb9 database, 294 completed (97.03%) BUSCOs were identified out of 303 groups searched, and 213 (70.30%) were single copy, 81 (26.73%) were duplicated, 7 (2.31%) were missing and only 2 (0.66%) fragmented (Additional file 4).


Coding regions were predicted for 46,576 sequences with GeneMarkS-T. Of the aligned sequences, 29.15% (8584) were classified as non-pine origin contigs, mainly belonging to the Fusarium genera (Fusarium mangiferae, F. fujikuroi, F. proliferatum, Fusarium nygamai, F. oxysporum). After filtering contaminants, non-frame selected and unannotated sequences, the final P. pinaster de novo transcriptome assembly contained 24,375 sequences, from which 17,593 (72.18%) were full-length genes. Best hit selection of BLAST alignment against the 4 databases generated a total of 20,864 (85.60%) unique sequences, from which 9957 (40.85%) were informative. EggNOG annotation predicted 23,674 sequences with family assignment, 5425 with at least one pathway (KEGG) assignment and 22,863 predicted protein domains. EggNOG associated 7614, 11,096 and 5050 gene ontology (GO) terms to biological process (BP), cellular compartment (CC) and molecular function (MF) categories, respectively. InterProScan predicted 20,188 protein domains and associated 6293, 1929 and 9346 GO terms to BP, CC and MF category, respectively (Additional file 5).

GhostKOALA assigned 12,741 K numbers (of 35,090) mainly classified according to the KEGG Orthology System in the functional categories genetic information processing, environmental information processing and cellular processes. Mercator assigned functional annotation to 15,347 P. pinaster transcripts (Additional file 6).

Mapping to the host and pathogen reference transcriptomes

Kallisto mapped a total of 964 million reads to the F. circinatum and P. pinaster combined dataset, which means 74.93% of the reads were mapped. In total, 72.97% of the reads were mapped to the P. pinaster de novo transcriptome reference, with similar mapping percentage between time points and between inoculated and mock-inoculated samples. Less than 8000 (< 0.01%) reads mapped to the F. circinatum reference transcriptome from any single mock-inoculated sample. For inoculated samples, 1.95% of the reads mapped to F. circinatum and an increase in the number of mapped reads was observed at the later time points owing to the growth of the pathogen (Additional file 7).

The principal component analysis (PCA) for P. pinaster and F. circinatum rlog data indicated clear separation of inoculated and mock-inoculated samples. The replicate samples show a high similarity with respect to the first two principal components for each time point. A small within group variance and a good separation of groups reflects the good quality of the analysis (Additional file 8).

Host and pathogen DE genes

A total of 13,323 differentially expressed (DE) genes were identified in inoculated P. pinaster samples compared to mock-inoculated ones. A notable increase of DE genes from 3 to 10 dpi was observed (Fig. 2). An unknown resistance protein gene was highly up-regulated at all time points and was the most up-regulated gene, with a log2 (Fold Change) of 25.08, 34.81 and 23.50 at 3, 5 and 10 dpi, respectively.

Fig. 2
figure 2

Venn’s diagram showing the overlap between differentially expressed (DE) genes. Left - Pinus pinaster DE genes at 3, 5 and 10 days post inoculation (dpi) in inoculated relative to mock-inoculated samples. Right – Fusarium circinatum DE genes between time points in inoculated samples. Red numbers – up-regulated genes, blue numbers – down-regulated genes, yellow numbers – genes contra-regulated between compared groups. 3v5: 3 relative to 5 dpi; 3v10: 3 relative to 10 dpi; 5v10: 5 relative to 10 dpi

For F. circinatum, at 3, 5 and 10 dpi 93.17% (4070 genes), 99.8% (4366 genes) and 99.9% (4372 genes) of the DE genes were considered high confident (HC) expressed genes, respectively (Additional file 7). When comparing DE genes at 3 versus 5 dpi in inoculated seedlings, 3427 genes were down-regulated and only 11 up-regulated. A similar pattern was observed between 3 versus 10 dpi with 4307 down-regulated and 21 up-regulated genes and between 5 and 10 dpi where 1599 genes were down-regulated and 24 up-regulated (Fig. 2; Additional file 7).

Over-represented GO terms in host DE gene clusters

Significant Pinus pinaster DE genes (13,323 genes) were classified in 53 clusters by using Hmisc R package. Due to the complexity, we set the |log2(Fold Change)| cut-off value to 1 (8802 DE genes) in order to reduce the number of clusters. DE genes were classified into 30 clusters, from which 12 had enriched GO terms and represent 92.7% of the DE genes (Fig. 3, Additional files 9 and 10).

Fig. 3
figure 3

Heatmap representing clusters with gene ontology (GO) enrichment for the significant Pinus pinaster differential expressed (DE) genes (FDR < 0.05; |log2(Fold Change inoculated / mock-inoculated)| = 1) at each time point (3, 5 and 10 days post-inoculation)

Genes in cluster 1 were up-regulated at 5 and 10 dpi and slightly at 3 dpi. Terms included in the BP category were related to response to stimulus, response to chitin, regulation of reactive oxygen species (ROS), response to oxidative stress, positive regulation of cell death and phosphorylation transduction, all responses commonly related to biotic stress. Phytohormone signaling was also evident, since terms related to ET, JA and SA were detected in cluster 1, as well as systemic acquire resistance (SAR) mediated by SA. In the MF category, terms were related to glycosyltransferase activity, catalytic activity, lipase activity and lyase activity. In the CC category, only two terms were enriched, plasma membrane and cell periphery.

Genes in cluster 2 were highly up-regulated at all time points with an increase from 3 to 10 dpi, and highly up-regulated at 3 dpi in comparison with the other clusters. In the BP category we found terms related to cellular metabolic processes, cellular localization, cellular component organization or biogenesis, cytosolic transport, protein-containing complex subunit organization, vesicle-mediated transport and proteasome-mediated ubiquitin-dependent protein catabolic process. MF terms were mostly related to binding activities, protein binding, heterocyclic compound binding, nucleotide binding and phosphorylase and hydrolase activities. CC terms were related to protein containing complex, vacuole membrane, membrane coat, vacuole coat and site of polarized growth.

Terms in cluster 3 in the BP category were related to flavonoid and anthocyanin metabolic processes, hormone transport, response to oxidative stress, response to JA and auxin polar transport. No terms were enriched in the MF category. NADH dehydrogenase complex, vacuole and mitochondria were found in the CC category. Genes of this cluster were mostly up-regulated at 10 dpi.

In cluster 5 only 3 BP terms were enriched, all related to transport. In cluster 11, a response to ethanol term was overrepresented in the BP category.

Genes classified in clusters 18, 19, 27 and 29 were down-regulated mostly at 10 dpi and slightly at 3 and 5 dpi (Fig. 3). Terms in these clusters were related to growth, development, reproduction, morphogenesis and photosynthesis. Terms related to isoprenoid, terpenoid and carotenoid metabolic processes were found.

Host phytohormone signaling

Four phytohormones seem to have a major role in P. pinaster defense response; they are JA, ET, SA and auxins. ABA, cytokinins (CK) and gibberellins seems to be suppressed. Log2(Fold Change) values for genes related to phytohormone signaling is provided as supplementary material (Additional file 11).

Jasmonic acid

Lipoxygenases (LOX) and 12-oxo-PDA-reductase genes (OPR) were up-regulated at all time points in inoculated seedlings. An allene oxidase cyclase gene (AOC) was only up-regulated at 10 dpi, while allene oxidase synthase genes (AOS) were up-regulated at 5 and 10 dpi but also down-regulated at 10 dpi. A receptor Coronatine Insensitive 1 (COI1) gene was down-regulated at 5 and 10 dpi. A JAZ (jasmonate ZIM domain) gene, a repressor of JA signaling, was highly up-regulated at all time points. A MYC transcription factor gene, which negatively regulates expression of JA responsive genes, was down-regulated at 5 and 10 dpi. TOPLESS (TPL) and NOVEL INTERACTOR of JAZ (NINJA) corepressors were down- and up-regulated at 10 dpi, respectively. Jasmonate methyl transferase (JMT) genes were up-regulated at all time points, indicating JA conversion to methyl jasmonate (MeJA) for systemic signaling. JA induced genes include PR proteins with roles in fungal cell wall degradation, such as chitinases (PR3) and β-1,3-glucanases (PR2) [27]. A total of 15 chitinase genes, mainly belonging to chitinase class VII, and a PR2 gene were up-regulated at all time points. Some peroxidases (PR9) can be induced by methyl jasmonate [26], and we detected up-regulation of 49 PR9 genes at all times, especially at 10 dpi. However, 14 PR9 genes were also down-regulated at 10 dpi (Fig. 4; Additional file 12).

Fig. 4
figure 4

Expression pattern of differential expressed (DE) genes related to pathogenesis related (PR) proteins in Pinus pinaster inoculated seedlings. For each PR protein, the average of the log2(Fold Change inoculated/ mock-inoculated) value is represented at 3, 5 and days post-inoculation. Error bars represent the standard error of the mean


In inoculated P. pinaster seedlings, 1-aminocyclepropane-1-carboxylic acid (ACC) synthase genes (ACS) were up-regulated at all time points while ACC-oxidase genes (ACO) were up-regulated at 5 and 10 dpi. The endoplasmic reticulum-associated receptor ETR2 was up-regulated at 5 and 10 dpi; however, ETR1 was down-regulated at 10 dpi. The downstream ETHYLENE INSENSITIVE 2 (EIN2), essential for positive regulation of ET signaling, was down-regulated at 5 and 10 dpi. EIN3 and EIL1 (EIN3-LIKE1), which act downstream of EIN2, were up-regulated at 5 and 10 dpi. ERF/AP2 ethylene-responsive transcription factors were up-regulated at all time points, especially at 5 and 10 dpi. EBF1 (EIN3-binding F-box protein 1), which degrades EIN3/EIL1 in the absence of ET, was down-regulated at 5 dpi.

Salicylic acid

Isochorismate synthase (ICS), the key enzyme involved in SA biosynthesis from the chorismate pathway, was down-regulated at 5 and 10 dpi. No differential expression of non-expressor of PR1 (NPR1) genes was found, however, genes of the TGA family transcription factors (TGA) were up-regulated at 10 dpi, suggesting SA signaling. SA can also be synthesized from the phenylalanine ammonia lyase (PAL) pathway, and PAL genes were up-regulated at all time points, although some PAL genes were also down-regulated. SA accumulation require two proteins, EDS1 (enhanced disease susceptibility 1) and PAD4 (Phytoalexin Deficient 4) [77]. Both EDS1 and PAD4 were up-regulated at 5 and 10 dpi. SA can be glycosylated by UDP-glycosyltransferase (UGT) to the inactive form 2–0-β-D-glucoside (SAG). Terms related to glycosyltransferase activity were enriched in cluster 1. PR1 is a PR protein commonly induced in defense response and its expression is SA responsive [107]. Three PR1 genes were up-regulated at all time points (Fig. 4). PR5 is SA, JA and ABA-responsive [113] and 10 PR5 genes were up-regulated at all time points, although two genes were also down-regulated (Fig. 4; Additional file 12).


Two YUCCA (indole-3-pyruvate monooxygenase) genes were down-regulated at all time points, although one YUCCA gene was also up-regulated at 10 dpi. TIR1 (Transport Inhibitor response), an auxin receptor [53], was up-regulated at 10 dpi. Two Aux/IAA (auxin/indole-3-acetic acid) family genes, which suppress the activity of transcriptional activators of the auxin response factor (ARF) family [99], were down-regulated at 10 dpi and 9 ARF genes were also down-regulated at 5 and 10 dpi. Small auxin up RNAs (SAURs), the largest family of auxin response genes, were up-regulated and down-regulated at 10 dpi, and several auxin-responsive genes were up- and down-regulated at all time points. GH3 genes, which inactivate IAA, were both up- and down-regulated at 5 and 10 dpi. The enzyme IAA carboxyl methyltransferase (IAMT) catalyzes the conversion of IAA into the inactive form methyl IAA (MeIAA) [78] and four IAMT1 genes were up-regulated at all dpi. Auxin efflux carriers PINFORMED (PIN) and P-glycoproteins (PGP) genes were down-regulated at 5 and 10 dpi, while an AUX3 gene, an auxin influx transporter related gene, was up-regulated at all time points.

Abscisic acid

Some 9-cis-epoxicarotenoid dioxygenase (NCED) genes, which participates in ABA biosynthesis, were up-regulated at all dpi, although one NCED gene was also down-regulated at 10 dpi. Zeaxanthin epoxidase (ZEP) was down-regulated at 10 dpi. Furthermore, ABA irreversible degradation was indicated since an ABA 8-hydroxylase gene was up-regulated at 5 and 10 dpi, which indicates ABA is converted to phaseic acid, with low activity. However, one ABA 8-hydroxylase gene was also down-regulated at 10 dpi. ABA signaling transduction is inhibited by type 2C protein phosphatases (PP2C), and PP2C genes were up-regulated at 5 and 10 dpi, although one gene was also down-regulated. Thus, ABA does not seem to have a major role in P. pinaster response to F. circinatum.

Gibberellic acid

Two key genes involved in GA biosynthesis, ent-kaurene synthase (KS) and ent-kaurene oxidase (KO), were down-regulated at 5 and 10 dpi. GA 20-oxidases (GA20ox) genes, involved in conversion of GA12 to the active forms of GA, was also down-regulated. GA 2-oxidases (GA2ox) convert bioactive gibberellins to their inactive form, and a GA2ox gene was highly up-regulated at all time points. These results suggest GA suppression.


We found down-regulation of histidine kinase receptor (HK) genes at 5 and 10 dpi. HK ultimately activates a family of transcription factors, ARR, though no changes in ARR genes were detected. A CK-O-glucosyltransferase gene was up-regulated at 5 and 10 dpi, which participate in CK degradation. Therefore, CK signaling seems to be suppressed in inoculated seedlings at 5 and 10 dpi, with no CK signaling activity at 3 dpi.


Some genes involved in BR biosynthesis were up-regulated at 5 and 10 dpi (DET2, BR60X) while others were down-regulated (DWF4 and DWF1). A receptor-like kinase BRI1 gene was down-regulated at 5 and 10 dpi. The downstream transcription factor BES1/BZR1 was up-regulated at 5 and 10 dpi. Brassinosteroid-responsive RING-H2 (BRH1) was also up-regulated at all time points, however, BRH1 is not only involved in response to BR stimulus, but also in response to chitin, and we found up-regulation of chitinase at all dpi.

Over-represented GO terms in pathogen DE gene clusters

Significant DE genes for F. circinatum (FDR < 0.05; |log2(FoldChange)| > 0.5) were classified into 7 clusters (Fig. 5; Additional file 9) and 2 of them had enriched GO terms (Additional file 13). These two clusters (cluster 3 and 4) represent 98% of the DE gene dataset and almost all of them (92.07%) were classified into cluster 3. The pattern of expression of both clusters was similar, increasing from 3 to 10 dpi (Fig. 5).

Fig. 5
figure 5

Heatmap representing clusters with gene ontology (GO) enrichment for Fusarium circinatum differentially expressed (DE) genes. Values are FPKM average at each time point (3, 5 and 10 days post-inoculation)

In cluster 3, enriched BP terms were related to generation of precursor metabolites and energy, such as mitochondrial ATP synthesis coupled electron transport, establishment of protein localization, transport, intracellular transport, oxidative phosphorylation and oxidation-reduction process. In the CC category enriched terms were related to mitochondrion, ribosome, cytoplasm, organelle membrane and endoplasmic reticulum. Phospholipid binding, phosphatidylinositol binding, lipid binding and structural constituent of ribosome were enriched in the MF category.

In cluster 4 enriched terms were related to perception and metabolic processes. In the BP category, enriched terms were related to regulation of cellular and metabolic processes, regulation of cyclin-dependent protein serine/threonine kinase activity, positive regulation of biological processes and signal transduction. Site of polarized growth was the only term enriched in the CC category. Osmosensor activity was enriched in the MF category.

Pathogen genes related to hormone production

Two genes related to GA biosynthesis (gibberellin cluster-C13-oxidase, gibberellin cluster-GA14-synthase) were up-regulated by F. circinatum at all time points. Genes that participate in geranyl geranyl diphosphate (GGDP) and ent-kaurene synthesis (gibberellin cluster-kaurensynthase and gibberellin cluster-GGPP-synthase) were also up-regulated, with an increase from 3 to 10 dpi (Fig. 6; Additional file 14).

Fig. 6
figure 6

Expression profile of Fusarium circinatum genes related to virulence. For each gene, FPKM values at each time point (3, 5 and 10 days post-inoculation) are indicated

Two genes (2-oxoglutarate-dependent ethylene succinate-forming enzyme and 2-keto-4-methylyhiobutyrate-dependent ethylene-forming enzyme) coding enzymes involved in fungal ET biosynthesis were up-regulated by F. circinatum at all time points, especially at the latter stage (Fig. 6). An auxin efflux gene was also up-regulated at all time points (Fig. 6, Additional file 14).

Isochorismatase hydrolase (ICSH) family proteins catalyze the hydrolysis of isochorismatase, a key metabolite for SA biosynthesis from the chorismate pathway. Two ICSH genes were up-regulated by F. circinatum and showed increased expression over time (Fig. 6, Additional file 14).

From the DE F. circinatum genes mentioned above, 4 of them showed hits to the Pathogen Host Interaction (PHI) database with E-value <1e− 4 (Additional file 15). Knockout of these genes in other pathogens resulted in reduced virulence in their hosts. Interestingly, one of these genes (2-keto-4-methylyhiobutyrate-dependent ethylene-forming enzyme) showed 90% similarity to a F. oxysporum gene, which when knocked out resulted in reduced virulence in tomato plants [74].


In forestry, sequencing of genomes is still a challenge, particularly for pines due to large genome sizes, retrotransposons and repeats. Therefore, transcriptome sequencing has become a good alternative for generating genomic resources, especially useful in understanding host-pathogen interactions [1, 114, 115]. In the present work, we generated a high quality de novo Pinus pinaster reference transcriptome constituted by 24,375 sequences, from which 17,523 were full length genes. Comparison of BUSCO results of the P. pinaster transcriptome with published Pinus transcriptomes showed high completeness and contiguity of the assembly [108, 110]. In total, 73% of reads mapped to the host transcriptome while 2% mapped to the pathogen, similar to results obtained in other studies [66]. Therefore, this high-quality shoot transcriptome represents a valuable resource for further research on maritime pine trees.

In addition, the present study provides a comprehensive characterization of the underlying molecular mechanism for defense and pathogenicity response in the Pinus pinaster - Fusarium circinatum pathosystem. The moderate resistance maritime pine has shown to the pathogen [49] can be explained, at least in part, by the early induction of defense-related genes and complex phytohormone signaling including SA, JA and ET. We also hypothesized key steps where the pathogen could be manipulating host defenses to its own benefit by altering host hormone homeostasis.

The early recognition and activation of P. pinaster defense responses can be deduced from genes classified in cluster 2, up-regulated at 3 dpi, before the fungus has penetrated within the host tissue. At 5 and 10 dpi, when the fungus has penetrated and invaded the host tissue, the expression of genes in cluster 1 highly increased and enriched GO terms related to phytohormone signaling, regulation of ROS, oxidative stress, positive regulation of cell death and signal transduction were found. This indicates activation of classical pattern triggered immunity (PTI) in response to the pathogen [11]. Terms related to chitinase activity were over-represented in this cluster and chitinases are commonly induced by pathogen attack in several trees such as P. abies, Pinus elliottii, Pinus sylvestris and Fagus sylvatica [27, 47, 71, 86, 88]. Rice and Arabidopsis perceive fungal chitin through the lysine motif (LysM) RLK CERK1 which induces CERK1 dimerization, essential for the activation of downstream signaling [67, 112]. We found up-regulation of genes encoding LysM motif RLK and chitinase in inoculated seedlings at 3 dpi, suggesting fungal recognition at the early stage. Furthermore, oligosaccharides released from chitin degradation can serve as pathogen-associated molecular patterns (PAMPs) for activation of PTI signaling in the host [68, 81].

Several studies have demonstrated the crucial role of phytohormones in host defense response [77]. In P. pinaster - F. circinatum pathosystem terms related to ET, JA and SA were over-represented in cluster 1. By contrast, GA, CK and ABA signaling seemed to be suppressed. Genes involved in JA and ET biosynthesis were up-regulated from 3 dpi and both phytohormones have been shown to cooperate under biotic stress conditions [9, 61, 62]. Despite the active biosynthesis of JA, the down-regulation of COI1 and MYC2 at 5 and 10 dpi, together with the up-regulation of JAZ genes reflect a block of JA signaling. COI1 has an F-box domain and is involved in the formation of Ubiquitin ligase E3 SCF complex, for protein degradation [30]. In the presence of JA-Ile, COI1 degrades JAZ proteins, which are direct targets of the SCF COI1 E3 ubiquitin-ligase [22, 98]. In the absence of COI1, JAZ proteins repress MYC transcription factor suppressing the expression of JA responsive genes. Based on the fact that COI1 is a key element for regulation of the JA signaling pathway [30] and the importance this phytohormone has shown in defense response against several necrotrophic pathogens [41], we hypothesize that F. circinatum could be blocking JA signaling by COI1 suppression. Similarly, Thatcher et al. [96] suggested that F. oxysporum hijacks COI1-mediated JA signaling to promote disease in Arabidopsis plants.

The induction of genes related to ET biosynthesis, together with the up-regulation of ethylene-responsive regulator genes ERF, ETR and EIN3/EIL1 suggests an active role of this phytohormone in P. pinaster defense response. However, the down-regulation of EIN2 at 5 and 10 dpi, essential for positive regulation of ET, could interfere with ethylene signaling. Interestingly, we found up-regulation of 2 genes related to ET biosynthesis in F. circinatum which could have a role in perturbing ET homeostasis in the host. Similarly, De Vleesschauwer et al. [28] reported that ET produced by the rice brown spot pathogen Cochliobolus miyabeanus is a virulence factor and suggested that fungal ET interferes with rice ET signaling to suppress effective defense pathways. Indeed, one of the two genes related to ET production in F. circinatum showed hits to the PHI database and knockout of this gene in F. oxysporum resulted in a reduced virulence phenotype in common bean plants [74].

ET and JA participate in induction of certain PR proteins [6, 26, 113]. In the present study, we detected up-regulation of PR1, PR2, PR3, PR5 and PR9 genes. Rakwal et al. [79] reported a coordinated increase in ET and some PR3 in rice plants. Similarly, Carrasco et al. [18] found synchronized increase between the induction of PR5 and ET in P. radiata seedlings inoculated with F. circinatum. PR2, PR3 and PR5 participate in degradation of glucans of the cell wall of the pathogen, making the fungus more susceptible to cell lysis and attack by other plant defense molecules [46, 107]. Similarly, during Fusarium culmorum infection, two basic isoforms of PR2 and three basic isoforms of PR3 were induced in germinating wheat seeds upon infection [19]. An increase of a class III PR9 responsible for lignin biosynthesis and cell wall thickening was detected in Pinus sylvestris roots infected with Heterobasidion annosum [1]. Indeed, induction of PR5 in P. radiata, highly susceptible to the pathogen, infected by F. circinatum, occurs after 6 dpi, but no differences were found at 3 dpi between inoculated and mock-inoculated samples [18]. Here we detected up-regulation since 3 dpi, before the pathogen has penetrated, suggesting quick activation of defense responses in P. pinaster against F. circinatum.

SA biosynthesis seems to occur from the PAL pathway and not from the chorismate pathway, since ICS was down-regulated at 10 dpi and not differentially expressed at 3 and 5 dpi. Similarly, Ding et al. [31] found induction of PAL in a wheat variety after Fusarium graminearum infection, while ICS was down-regulated, suggesting SA accumulation via the phenylpropanoid pathway. Interestingly, we found up-regulation of two ICSH genes by F. circinatum, which catalyzes the hydrolysis of ICS. Moreover, ICSH1 has been shown to accumulate in a highly aggressive Verticillium dahliae isolate but not in a weakly aggressive isolate [34]. Indeed, cotton plants inoculated with a V. dahliae mutant lacking the ICSH1 gene resulted in attenuated aggressiveness while higher levels of SA were detected [57]. Zhu et al. [123] propose that ICSH1 is a virulence factor in V. dahliae and that the high level of ICSH in highly aggressive isolates represses the SA pathway in potato plants by hydrolyzing ICS. Similarly, we suggest that F. circinatum prevents SA biosynthesis from the chorismate pathway, although whether ICSH is a virulence factor in F. circinatum requires further investigation. In addition, TGA, EDS1 and PAD4 were up-regulated in inoculated P. pinaster seedlings, suggesting SA accumulation has a role in the host defense response. Furthermore, PAD4 is also required for production of the phytoalexin camalexin and PR1 synthesis. PR1 gene expression is SA responsive and three PR1 genes were up-regulated at all dpi. Carrasco et al. [18] also found PR1 genes highly up-regulated in F. circinatum resistant P. radiata genotypes.

Van der Does et al. [106] proposed that suppression of the JA pathway by SA functions downstream of the E3 ubiquitin-ligase Skip-Cullin-F-box complex SCFCOI1. However, cooperation between JA and SA has been also reported [36, 95]. Due to the up-regulation of genes related to SA and JA biosynthesis in inoculated P. pinaster seedlings and based on the fundamental role JA has shown against necrotrophic pathogens [41], we suggest cooperation of SA and JA rather than an antagonism in the P. pinaster-F. circinatum interaction.

Genes encoding enzymes that catalyze the methylation of JA (MeJA) were up-regulated at all time points suggesting a role of systemic signaling in P. pinaster defense response. Furthermore, “SAR mediated by SA” term was over-represented in cluster 1. Truman et al. [102] suggested that JA has an earlier role in SAR establishment, before SA accumulation in Arabidopsis. Similarly, we found up-regulation of genes related to JA biosynthesis (LOX, OPR) since 3 dpi, before induction of SA related genes, mostly up-regulated at 5 and 10 dpi (NPR1, TAG, PAD4, EDS1), suggesting again, possible cooperation between the two phytohormones. Foliar application of MeJA in P. radiata conferred increased resistance to the necrotroph Diplodia pinea [42] and against other pathogenic fungi in P. abies [39, 54]. However, MeJA pre-treatment was ineffective in protection against F. circinatum in Pinus patula [37] and P. pinaster seedlings [111]. Nevertheless, the high induction of MeJA producing genes in P. pinaster inoculated seedlings indicates a role in defense response against the pathogen. Indeed, Sasaki et al. [85] described MeJA as a key component in the jasmonate signaling pathway by controlling its own expression via a feedback mechanism by inducing the expression of LOX and AOS, both enzymes that catalyze key steps in JA biosynthesis. This can explain the up-regulation of JA biosynthesis genes in spite of the suppression of JA signaling. Furthermore, studies have reported the role of MeJA in chalcone synthase (CHS) induction in soybean and parsley [25] and Picea glauca [82]. CHS is a key enzyme in the flavonoid biosynthesis pathway and flavonoids play an important role in plant defense against pathogens. We found up-regulation of CHS at all dpi in P. pinaster inoculated seedlings, with an increase from 3 to 10 dpi.

Although SA, JA and ET are key players in P. pinaster response to F. circinatum, auxins can also have a role since some SAUR genes and a TIR1 gene were up-regulated. However, we also detected suppression of auxin biosynthesis and signaling (down-regulation of indole-3-pyruvate monooxygenase and ARF), as well as accumulation and conjugation to inactive forms inside the cell (up-regulation of GH3 and a AUX gene encoding for influx proteins, and down-regulation of PIN and PGP genes encoding efflux proteins). Nevertheless, IAMT1 genes were up-regulated suggesting conversion of IAA into the inactive MeIAA form. Crosstalk between auxins and other phytohormones has been documented and an inhibitory effect of auxins on JA signaling in Arabidopsis seedlings has been reported [56]. Moreover, auxins are inducers of expansins, involved in cell wall extensibility [20] which could enhance pathogen penetration. Indeed, auxins mostly are related to growth and development which are functions sacrificed by the plant to favor mechanisms involved in defense responses. This trade-off between growth and defense is reflected by over-represented GO terms related to growth, morphogenesis and photosynthesis for clusters with down-regulated DE genes (clusters 18, 19, 27 and 29; Fig. 3). Nonetheless, the role of auxins in P. pinaster defense response to F. circinatum is not clear and needs to be elucidated.

Furthermore, a gene related to auxin efflux was up-regulated by F. circinatum at all time points, suggesting auxin accumulation could have a role in virulence. Similarly, an over-accumulation of IAA in F. oxysporum was related to a hypervirulent phenotype on Orobanche [24]. The two-fold reduction of a gene required for auxin biosynthesis in Puccinia graminis f. sp. tritici led to a decrease in pustule formation [119]. GA production is common in species of the FFC. Although Malonek et al., [64] reported that only one gene in the GA biosynthetic cluster is present in F. circinatum, we found two potential genes, related to GA biosynthesis, up-regulated. These two genes also had hits to the PHI database, and knockout of these genes in other pathogens resulted in reduced virulence and loss of pathogenicity phenotypes in the host. However, these results must be considered with caution, since percentages of identity ranged from 22 to 29%.

By 10 days after inoculation, when the fungus has invaded the host tissue and the lesion in the shoot tip is evident, an additional defense response seems to occur in P. pinaster seedlings (cluster 3). GO terms related to flavonoids and anthocyanins biosynthesis, secondary metabolites involved in defense responses, were over-represented [7, 100]. Despite the effort of the plant in synthesizing these secondary metabolites at the latter stage of infection, the host is not able to counteract fungal infection.

This study highlights the importance of phytohormones in the P. pinaster-F. circinatum interaction. One possibility would be to use hormone application to induce resistance. While previous results show the application of MeJA at various concentrations does not affect lesion development of F. circinatum inoculated P. pinaster seedlings [111], this study points to SA in combination with MeJA as a potential strategy to investigate.


This work provides knowledge of mechanisms underlying the P. pinaster defense response against F. circinatum, indicating activation of defense mechanisms from as early as 3 dpi, by induction of PR genes and mainly regulated by a complex signaling pathway involving crosstalk between SA, JA and ET. Moreover, we hypothesize key steps where the pathogen could be manipulating host defense in its favor, mainly by perturbing phytohormone signaling homeostasis in the host. Future work in measuring SA, JA, ET and auxin content in planta, as well as functional studies with F. circinatum mutants, will be necessary to support this hypothesis.


Plant and fungal material

Six-month-old P. pinaster seedlings purchased at ‘Eskalmendi’ nursery (Alava, Spain) were used for the experiment. These seedlings were grown in this nursery from seeds of P. pinaster, provenance Landas. Seedlings were maintained in a greenhouse at 20–22 °C with a photoperiod of 12 h light / 12 h darkness and inoculated after 2 weeks of acclimation.

For inoculations, a virulent F. circinatum isolate from Basque Country (Northern Spain) (Isolate CECT20759, isolated from a P. radiata tree) [50] was used. A fungal spore suspension in sterile distilled water was prepared after 1 week of culture on Potato Dextrose Agar (PDA) by scraping the plate surface and passing through two layers of glass wool. Spore concentration was measured with a hemocytometer and adjusted to 5 × 105 spores/ml.

Inoculation and microscopic observation

Ninety six-month-old P. pinaster seedlings were inoculated with a green fluorescent protein (GFP)-tagged strain of F. circinatum. Agrobacterium-mediated transformation of isolate CECT20759-GFP was performed as previously described [69, 80] using the pRF-gGFP plasmid [75]. This plasmid contained the sGFP coding sequence under the control of the A. nidulans Pgpd promoter excised from plasmid gGFP [65]. The first 2 cm of the shoot tip were excised and a 2 μl drop of the spore suspension (1000 conidia) was deposited in the wound with a micropipette. A set of 54 seedlings were mock-inoculated with sterile distilled water. Plants were covered with a plastic bag for 24 h in order to maintain high humidity and favor fungal infection. During the first 8 days following inoculation, as well as at 10, 14 and 17 dpi, plant tissue from 6 inoculated and 2 mock-inoculated plants per day were visualized under an epifluorescence microscope (Nikon EFD-3). Cross-sections were placed on glass slides, covered with cover glasses and visualized under the epifluorescence microscope. Sections were taken from the point of inoculation and progressively downward until no pathogen was visualized, a minimum of four sections were visualized for each plant. Progression of pathogen growth within host tissue served to determine the times of sampling for RNAseq analyses.

Inoculation and tissue sampling for RNA extraction

A total of 288 P. pinaster seedlings were used for the experiment. Half of them were inoculated with the pathogen as explained above, while the remaining half were mock-inoculated with sterile distilled water. For sampling, the top 1.5 cm of shoot tissue was harvested for each seedling at three different times: 3, 5, and 10 dpi, for both inoculated and mock-inoculated seedlings. We used 4 biological replicates (BR) per group (mock-inoculated and inoculated) with a pool of 8 individuals each. Plant material was immediately frozen in liquid nitrogen and stored at − 80 °C until use.

A total of 16 inoculated and 16 non-inoculated seedlings were maintained for visualizing disease progress. Lesion length was measured at the end of the experiment (33 dpi). Additionally, 3 cm of the tip of 6 randomly collected plants of each group (inoculated and mock-inoculated) were surface sterilized by immersion in 70% EtOH for 1 min. Plant tissue was transversally cut and cultured on Fusarium Selective Medium [2] to verify efficacy of the inoculation.

RNA isolation and sequencing

Total RNA was extracted using a Plant/Fungi Total RNA Purification Kit (Norgen Biotek Corp., Thorold, Ontario) following the manufacturer’s instructions and stored at − 80 °C. The protocol included a DNase treatment step for removal of residual DNA (Norgen’s RNase-free DNase I Kit). The integrity of extracted RNA was assessed using a 2100 Bioanalyzer (Agilent Technologies). RNA with an RNA Integrity Number (RIN) > 7 was considered good quality. Two biological replicates for inoculated samples at 10 dpi did not pass the RIN threshold and were excluded for the rest of the analysis, likely due to necrosis of the tissue (Additional file 16).

Approximately 40 μl of total RNA for each sample was submitted to Macrogen (Macrogen Korea). Twenty-four TruSeq mRNA stranded libraries with polyA selection and an insert size of 300 bp was prepared. One hundred one bp paired-end reads were generated with Illumina HiSeq 4000, with a sequencing depth of 81–181 million reads per sample. Samples were multiplexed in 4 lanes, with 6 samples per lane. In order to avoid technical errors due to sample position in the sequencer, each biological replicate of each sample was included in a different lane.

Raw data quality control and filtering

Quality control of raw Ilumina reads was performed by FastQC analyses v0.11.7 [3]. Trimmomatic 0.36 was used for trimming and filtering of low-quality reads and adapter removal [10]. Bases with a Phred Score below 30 and reads shorter than 40 bp were removed. The first 15 bases of all reads were trimmed to remove sequencing biases. The quality of trimmed and filtered reads was checked again by FastQC (Additional file 16).

Reference transcriptomes

Pinus pinaster de novo transcriptome assembly

Preliminary assemblies

Diverse studies have shown a high variability in the assembly of reads when using different assemblers [13, 122]. Therefore, the de novo transcriptome assembly of P. pinaster was created using two different assemblers: Trinity v.2.4.0 [43] and transABySS v.2.0.1 [83], programs specially developed for de novo transcriptome assembly of RNAseq short-read data. A multi kmer strategy was adopted; De Bruijn graphs were built over different kmer values and the resulting assemblies merged to improve sensitivity and accuracy of gene set reconstruction [33, 93, 122].

Seven assemblies were created with Trinity with kmer values of 19, 21, 23, 25, 27, 29 and 31 and the next set of parameters: minimum contig length of 350 bp and maximum read coverage of 50 for in silico normalization. A further 5 assemblies were generated with kmer sizes of 19, 21, 23, 25 and 27 without normalization of the reads. For TransABySS, we ran 9 different de novo assemblies with a kmer value ranging from 21 to 77, with a step size of 8, as well as a kmer value of 25. The minimum output sequence length was set to 350 bp for all TransABySS assemblies. Quality of each preliminary assembly was checked using Transrate v.1.0.3 [91].

Merging assemblies

Best quality preliminary assemblies were merged into one dataset using the EvidentialGene tr2aacds pipeline version 2017.12.21 [40], which builds the optimal assembly from a pool of different assemblies. It uses fastanrdb of exonerate package v 2.2.0 [90] for pairwise sequence comparison based on protein qualities for predicting the best coding DNA sequences (CDS) among identical sequences and reduce redundancy. Then, cd-hit-est v.4.7 [55] and BLASTn (blast v.2.7.0) [121] are used to cluster nucleotide sequences with 98% similarity into loci. EvidentialGene tr2aacds returns the subset of most accurate coding genes, classified into alternate or main (pimary) CDS. Alternate CDS for each gene were discarded. Additionally, an “okay” or “drop” value is assigned using scores of alignment and protein quality to separate “useful” and “not useful” transcripts, respectively. The “drop” class contains redundant and uninformative mRNA transcripts and they were discarded for the rest of the analyses. The quality of the transcript sequences derived from the merged assembly were checked by Transrate v.1.0.3.


Both reference transcriptomes were annotated using the Eukaryotic Non-Model Transcriptome Annotation Pipeline (EnTap) version 0.8.2 [44]. Coding regions of the transcripts were selected by GenemarkS-T v5.1 March 2014 [94]. The annotation process integrates similarity search across different databases with a minimum query and target coverage of 80 and 60%, respectively. NCBI non-redundant protein (release 2018–03), RefSeq (release 87), SwissProt (release 2018–03) and Arabidopsis proteome (release 2018.03) databases were used for BLASTp alignment using Diamond 0.9.9 [16]. Non-pine origin sequences were removed from the assembly by significant alignment to fungal, bacterial, viral, insect, archaea, opistokhonta and amoebozoa sequences. For orthologous group and gene ontology (GO) assignment InterProScan v5.28–67.0 [51] and EggNOG v0.12.7 [48] were used. Finally, all contaminants and non-frame selected and unannotated sequences were manually filtered. Kyoto Encyclopaedia of Genes and Genomes (KEGG) orthology (KO) annotation was also assigned by using GhostKOALA [52].

Mercator [60], with default parameters and including all the databases available, was used to assign predicted proteins into MapMan v.3.5.1R2 [97] bins, a tool that allows the visualization of metabolic pathways and processes of a large set of data.

Assembly validation

Completeness and contiguity of the assembly was checked using BUSCO v3.0.2 (Benchmarking of Single-Copy Orthologs) [89]. We used the eukaryote_odb9 and embryophyta_odb9 lineages to identify putative universal single copy orthologs (USCOs) in the assembly.

Fusarium circinatum reference transcriptome

The F. circinatum reference transcriptome was obtained from the F. circinatum (strain FSP34) genome sequence (B. D [117].) by extracting the longest transcript sequence for all predicted genes (15,049). A total of 14,185 sequences were annotated with EnTAP, of which 5368 were assigned GO terms [109]

Mapping, differential expression (DE) and gene ontology (GO) enrichment analysis

The F. circinatum transcriptome and P. pinaster de novo transcriptome assembly were combined in a single dataset to account for cross species mapping [72]. This dataset was used to map reads with Kallisto v.0.44.0 [15] with sequence bias correction and bootstrap samples set to 100. Kallisto abundance output files (transcript abundance estimates) for each read was imported to R 3.5.1 with tximport v.1.6.0 [92]. DESeq2 v.1.18.1 [63] was used for DE analysis, with a FDR < 0.05 and |log2(Fold change)| > 0.5. Transcripts with less than 20 reads in at least 3 samples were filtered for the rest of the analysis. For host DE analysis F. circinatum data was removed and expression levels at different time points (3, 5 and 10 dpi) were compared between inoculated and mock-inoculated samples. For F. circinatum DE analyses, P. pinaster data was filtered and gene expression levels between time points were compared between the inoculated samples (3 versus 5 dpi, 3 versus 10 dpi and 5 versus 10 dpi). To excluded potential endophyte contamination and confirm the confidence of F. circinatum expressed genes, a high confidence DE analysis (inoculated relative to mock-inoculated) with all the data (host and pathogen) was performed and then compared to the one with only inoculated samples. Genes considered not high confident at all time points were discarded. GOSeq v.1.34.0 [120] was used to identify GO terms significantly overrepresented (FDR < 0.10) in the DE data. GO enrichment was based on the annotated transcriptomes for each species.

PCA was performed for P. pinaster and F. circinatum normalized read counts (regularized-logarithm transformation or rlog - DESeq2 rlog function) for all samples, to visualize the overall effect of experimental covariates and batch effects.

Significant DE genes for each species were clustered based on the FPKM (Fragments Per Kilobase of transcript per Million mapped reads) values using Hmisc v.4.1–1, with a Pearson correlation of 0.5. Overrepresented GO terms for each cluster were identified with GOSeq. F. circinatum genes were subjected to a PHI-BLAST analysis using the PHI database 4.2 [104] to identify potential pathogenicity and virulence factors. Default parameters were employed and lowest e-value hits were considered further.

Availability of data and materials

The datasets generated and analysed during the current study are available in the Sequence Read Archive (SRA) repository, accessible through BioProject accession PRJNA543723.



Abscisic acid


1-aminocyclepropane-1-carboxylic acid




ACC synthase


Allene oxidase cyclase


Allene oxidase synthase


Auxin response factor


Auxin/indole-3-acetic acid


Biological process


Biological replicate


Cellular compartment


Chalcone synthase




Coronatine Insensitive 1


Differentially expressed


Days post-inoculation


EIN3-binding F-box protein 1


Enhanced disease susceptibility 1


Ethylene insensitive 2




Endoplasmic reticulum-associated receptors


Fragments Per Kilobase of transcript per Million mapped reads


Gibberellic acid


GA 20-oxidases


GA 2-oxidases


Green fluorescent protein


Geranyl geranyl diphosphate


Gretchen Hagen 3


Gene ontology


High confident


Histidine kinase


Indol-3-acetic acid




IAA carboxyl methyltransferase


Isochorismate synthase


Isochorismatase hydrolase


Jasmonic acid




Jasmonate ZIM domain


Jasmonate methyl transferase


Ent-kaurene oxidase


Ent-kaurene synthase




Methyl IAA


Methyl jasmonate


Molecular function


9-cis-epoxicarotenoid dioxygenase


Novel interactor of JAZ


Non-expressor of PR1




Phytoalexin deficient 4


phenylalanine ammonia lyase


Pathogen-associated molecular patterns


Principal component analysis


Phosphorylase family protein




pathogen-host interactions




2C protein phosphatases




Pattern triggered immunity


RNA Integrity Number


Reactive oxygen species


Salicylic acid






Systemic acquire resistance


Small auxin up RNAs


Transport Inhibitor response






Zeaxanthin epoxidase


  1. Adomas A, Heller G, Li G, Olson Å, Chu TM, Osborne J, Craig D, et al. Transcript profiling of a conifer pathosystem: response of Pinus sylvestris root tissues to pathogen (Heterobasidion annosum) invasion. Tree Physiol. 2007;27(10):1441–58.

    Article  CAS  PubMed  Google Scholar 

  2. Aegerter BJ, Gordon TR. Rates of pitch canker induced seedling mortality among Pinus radiata families varying in levels of genetic resistance to Gibberella circinata (anamorph Fusarium circinatum). For Ecol Manag. 2006;235(1–3):14–7.

    Article  Google Scholar 

  3. Andrews, S. (2012) FastQC a quality control tool for high throughput sequence data. Retrieved from

    Google Scholar 

  4. Arnerup J. Induced defence responses in Picea abies triggered by Heterobasidion annosum s . l. Uppsala: Swedish University of Agricultural Sciences; 2011. 56 p

    Google Scholar 

  5. Bari R, Jones JDG. Role of plant hormones in plant defence responses. Plant Mol Biol. 2009;69(4):473–88.

    Article  CAS  PubMed  Google Scholar 

  6. Belhadj A, Telef N, Saigne C, Cluzet S, Barrieu F, Hamdi S, Mérillon JM. Effect of methyl jasmonate in combination with carbohydrates on gene expression of PR proteins, stilbene and anthocyanin accumulation in grapevine cell cultures. Plant Physiol Biochem. 2008;46(4):493–9.

    Article  CAS  PubMed  Google Scholar 

  7. Bennett RN, Wallsgrove RM. Secondary metabolites in plant defence mechanisms. New Phytol. 1994;127(4):617–33.

    Article  CAS  PubMed  Google Scholar 

  8. Berrocal-Lobo M, Molina A. Ethylene response factor 1 mediates Arabidopsis resistance to the Soilborne fungus Fusarium oxysporum. Mol Plant-Microbe Interact. 2004;17(7):763–70.

    Article  CAS  PubMed  Google Scholar 

  9. Berrocal-Lobo M, Molina A, Solano R. Constitutive expression of ETHYLENE-RESPONSE-FACTOR1 in Arabidopsis confers resistance to several necrotrophic fungi. Plant J. 2002;23(1):23–32.

    Article  Google Scholar 

  10. Bolger AM, Lohse M, Usadel B. Genome analysis trimmomatic : a flexible trimmer for Illumina sequence data. Bioinformatics. 2014;30(15):2114–20.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  11. Boller T, He SY. Innate immunity in plants : an arms race between pattern recognition receptors in plants and effectors in microbial pathogens. Science. 2009;324(5928):742–4.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  12. Bömke C, Tudzynski B. Diversity, regulation, and evolution of the gibberellin biosynthetic pathway in fungi compared to plants and bacteria. Phytochemistry. 2009;70(15–16):1876–93. Elsevier Ltd.

    Article  CAS  PubMed  Google Scholar 

  13. Bradnam KR, Fass JN, Alexandrov A, Baranay P, Bechner M, Birol I, Boisvert S, et al. Assemblathon 2 : evaluating de novo methods of genome assembly in three vertebrate species. Gigascience. 2013;2(10):1–31.

    Google Scholar 

  14. Bragança H, Diogo E, Moniz F, Amaro P. First report of pitch canker on pines caused by Fusarium circinatum in Portugal. Plant Dis. 2009;93(10):1079.

    Article  PubMed  Google Scholar 

  15. Bray NL, Pimentel H, Melsted P, Pachter L. Near-optimal probabilistic rna-seq quantification. Nat Biotechnol. 2016;34(5):4–8.

    Article  CAS  Google Scholar 

  16. Buchfink B, Xie C, Huson DH. Fast and sensitive protein alignment using DIAMOND. Nat Methods. 2014;12(1):59–60.

    Article  CAS  PubMed  Google Scholar 

  17. Canales J, Bautista R, Label P, Gómez-Maldonado J, Lesur I, Fernández-Pozo N, Rueda-López M, et al. De novo assembly of maritime pine transcriptome: implications for forest breeding and biotechnology. Plant Biotechnol J. 2014;12(3):286–99.

    Article  CAS  PubMed  Google Scholar 

  18. Carrasco A, Wegrzyn JL, Durán R, Fernández M, Donoso A, Rodriguez V, Neale D, et al. Expression profiling in Pinus radiata infected with Fusarium circinatum. Tree Genet Genomes. 2017;13(2):46.

    Article  Google Scholar 

  19. Caruso C, Chilosi G, Caporale C, Leonardi L, Bertini L, Magro P, Buonocore V. Induction of pathogenesis-related proteins in germinating wheat seeds infected with Fusarium culmorum. Plant Sci. 1999;140(1):87–97.

    Article  CAS  Google Scholar 

  20. Catalá C, Rose JKC, Bennett AB. Auxin-regulated genes encoding cell wall-modifying proteins are expressed during early tomato fruit growth. Plant Physiol. 2000;122(2):527–34.

    Article  PubMed  PubMed Central  Google Scholar 

  21. Chanclud E, Morel JB. Plant hormones: a fungal point of view. Mol Plant Pathol. 2016;17(8):1289–97.

    Article  PubMed  PubMed Central  Google Scholar 

  22. Chini A, Fonseca S, Fernández G, Adie B, Chico JM, Lorenzo O, García-Casado G, et al. The JAZ family of repressors is the missing link in jasmonate signalling. Nature. 2007;448(7154):666–71.

    Article  CAS  PubMed  Google Scholar 

  23. Clarke JD, Volko SM, Ledford H, Ausubel FM, Dong X. Roles of salicylic acid, jasmonic acid, and ethylene in cpr-induced resistance in Arabidopsis. Plant Cell. 2000;12(11):2175–90.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  24. Cohen BA, Amsellem Z, Maor R, Sharon A, Gressel J. Transgenically enhanced expression of indole-3-acetic acid confers hypervirulence to plant pathogens. Phytopathology. 2002;92(6):590–6.

    Article  CAS  PubMed  Google Scholar 

  25. Creelman RA, Tierneyt ML, Mullet JE. Jasmonic acid/methyl jasmonate accumulate in wounded soybean hypocotyls and modulate wound gene expression (chalcone synthase/vegetative storage protein/prolne-rich cell wail protein). Plant Biol. 1992;89:4938–41 Retrieved from

    CAS  Google Scholar 

  26. Curtis MD, Rae AL, Rusu AG, Harrison SJ, Manners JM. A peroxidase gene promoter induced by phytopathogens and methyl jasmonate in transgenic plants. Mol Plant-Microbe Interact. 1997;10(3):326–38.

    Article  CAS  PubMed  Google Scholar 

  27. Davis JM, Wu H, Cooke JEK, Reed JM, Luce KS, Michler CH. Pathogen challenge , salicylic acid , and jasmonic acid regulate expression of chitinase gene homologs in pine. Mol Plant-Microbe Interact. 2002;15(4):380–7.

    Article  CAS  PubMed  Google Scholar 

  28. De Vleesschauwer D, Yang Y, Vera Cruz C, Hofte M. Abscisic acid-induced resistance against the brown spot pathogen Cochliobolus miyabeanus in rice involves MAP kinase-mediated repression of ethylene signaling. Plant Physiol. 2010;152(4):2036–52.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  29. Desjardins AE, Manandhar HK, Plattner RD, Manandhar GG, Poling SM, Maragos CM. Fusarium species from nepalese rice and production of mycotoxins and gibberellic acid by selected species. Appl Environ Microbiol. 2000;66(3):1020–5.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  30. Devoto A, Nieto-Rostro M, Xie D, Ellis C, Harmston R, Patrick E, Davis J, et al. COI1 links jasmonate signalling and fertility to the SCF ubiquitin-ligase complex in Arabidopsis. Plant J. 2002;32(4):457–66.

    Article  CAS  PubMed  Google Scholar 

  31. Ding L, Xu H, Yi H, Yang L, Kong Z, Zhang L, Xue S, et al. Resistance to hemi-biotrophic F graminearum infection is associated with coordinated and ordered expression of diverse defense signaling pathways. PLoS One. 2011;6(4):e19008.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  32. Donoso A, Rodriguez V, Carrasco A, Ahumada R, Sanfuentes E, Valenzuela S. Relative expression of seven candidate genes for pathogen resistance on Pinus radiata infected with Fusarium circinatum. Physiol Mol Plant Pathol. 2015;92:42–50. Elsevier Ltd.

    Article  CAS  Google Scholar 

  33. Durai DA, Schulz MH. Informed kmer selection for de novo transcriptome assembly. Bioinformatics. 2016;32(11):1670–7.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  34. El-Bebany AF, Rampitsch C, Daayf F. Proteomic analysis of the phytopathogenic soilborne fungus Verticillium dahliae reveals differential protein expression in isolates that differ in aggressiveness. Proteomics. 2010;10(2):289–303.

    Article  CAS  PubMed  Google Scholar 

  35. Elvira-Recuenco M, Iturritxa E, Majada J, Alia R, Raposo R. Adaptive potential of maritime pine (Pinus pinaster) populations to the emerging pitch canker pathogen, Fusarium circinatum. PLoS One. 2014;9(12):1–21.

    Article  CAS  Google Scholar 

  36. Engelberth J, Viswanathan S, Engelberth MJ. Low concentrations of salicylic acid stimulate insect elicitor responses in Zea mays seedlings. J Chem Ecol. 2011;37(3):263–6.

    Article  CAS  PubMed  Google Scholar 

  37. Fitza KNE, Payn KG, Steenkamp ET, Myburg AA, Naidoo S. Chitosan application improves resistance to Fusarium circinatum in Pinus patula. S Afr J Bot. 2013;85:70–8. South African Association of Botanists.

    Article  CAS  Google Scholar 

  38. Flaishman MA. Timing of fungal invasion using host’s ripening hormone as a signal. Proc Natl Acad Sci. 1994;91(14):6579–83.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  39. Franceschi VR, Krekling T, Christiansen E. Application of methyl jasmonate on Picea abies (Pinaceae) stems induces defense-related responses in phloem and xylem. Am J Bot. 2002;89(4):578–86.

    Article  CAS  PubMed  Google Scholar 

  40. Gilbert, D. (2013) EvidentialGene: tr2aacds, mRNA transcript assembly software. Retrieved from

    Google Scholar 

  41. Glazebrook J. Contrasting mechanisms of defense against biotrophic and necrotrophic pathogens. Annu Rev Phytopathol. 2005;43(1):205–27.

    Article  CAS  PubMed  Google Scholar 

  42. Gould N, Reglinski T, Spiers M, Taylor JT. Physiological trade-offs associated with methyl jasmonate - induced resistance in Pinus radiata. Can J For Res. 2008;38(4):677–84.

    Article  Google Scholar 

  43. Grabherr MG, Haas BJ, Yassour M, Levin JZ, Thompson DA, Amit I, Adiconis X, et al. Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nat Biotechnol. 2011;29(7):644–52.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  44. Hart AJ, Ginzburg S, Xu MS, Fisher CR, Rahmatpour N, Mitton JB, Paul R, et al. EnTAP: bringing faster and smarter functional annotation to non-model eukaryotic transcriptomes. bioRxiv. 2018:307868.

  45. Hayden KJ, Garbelotto M, Knaus BJ, Cronn RC, Rai H, Wright JW. Dual RNA-seq of the plant pathogen Phytophthora ramorum and its tanoak host. Tree Genet Genomes. 2014;10(3):489–502.

    Article  Google Scholar 

  46. Hématy K, Cherk C, Somerville S. Host-pathogen warfare at the plant cell wall. Curr Opin Plant Biol. 2009;12(4):406–13.

    Article  CAS  PubMed  Google Scholar 

  47. Hietala AM, Kvaalen H, Schmidt A, Jøhnk N, Solheim H, Fossdal CG. Temporal and spatial profiles of chitinase expression by Norway spruce in response to bark colonization by Heterobasidion annosum. Appl Environ Microbiol. 2004;70(7):3948–53.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  48. Huerta-Cepas J, Szklarczyk D, Forslund K, Cook H, Heller D, Walter MC, Rattei T, et al. EGGNOG 4.5: a hierarchical orthology framework with improved functional annotations for eukaryotic, prokaryotic and viral sequences. Nucleic Acids Res. 2016;44(D1):D286–93.

    Article  CAS  PubMed  Google Scholar 

  49. Iturritxa E, Ganley RJ, Raposo R, García-Serna I, Mesanza N, Kirkpatrick SC, Gordon TR. Resistance levels of Spanish conifers against Fusarium circinatum and Diplodia pinea. For Pathol. 2013;43(6):488–95.

    Article  Google Scholar 

  50. Iturritxa E, Ganley RJ, Wright J, Heppe E, Steenkamp ET, Gordon TR, Wingfield MJ. A genetically homogenous population of Fusarium circinatum causes pitch canker of Pinus radiata in the Basque Country, Spain. Fungal Biol. 2011;115(3):288–95.

    Article  PubMed  Google Scholar 

  51. Jones P, Binns D, Chang HY, Fraser M, Li W, McAnulla C, McWilliam H, et al. InterProScan 5: genome-scale protein function classification. Bioinformatics. 2014;30(9):1236–40.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  52. Kanehisa M, Sato Y, Morishima K. BlastKOALA and GhostKOALA: KEGG tools for functional characterization of genome and metagenome sequences. J Mol Biol. 2016;428(4):726–31. The authors.

    Article  CAS  PubMed  Google Scholar 

  53. Kepinski S, Leyser O. The Arabidopsis F-box protein TIR1 is an auxin receptor. Nature. 2005;435(7041):446–51.

    Article  CAS  PubMed  Google Scholar 

  54. Kozlowski G, Buchala A, Métraux JP. Methyl jasmonate protects Norway spruce [Picea abies (L.) karst.] seedlings against Pythium ultimum Trow. Physiol Mol Plant Pathol. 1999;55(1):53–8.

    Article  CAS  Google Scholar 

  55. Li W, Godzik A. Cd-hit: a fast program for clustering and comparing large sets of protein or nucleotide sequences. Bioinformatics. 2006;22(13):1658–9.

    Article  CAS  PubMed  Google Scholar 

  56. Liu J, Wang XJ. An integrative analysis of the effects of auxin on jasmonic acid biosynthesis in Arabidopsis thaliana. J Integr Plant Biol. 2006;48(1):99–103.

    Article  CAS  Google Scholar 

  57. Liu T, Song T, Zhang X, Yuan H, Su L, Li W, Xu J, et al. Unconventionally secreted effectors of two filamentous pathogens target plant salicylate biosynthesis. Nat Commun. 2014;5:4686.

    Article  CAS  PubMed  Google Scholar 

  58. Llorente F, Muskett P, Sánchez-Vallet A, López G, Ramos B, Sánchez-Rodríguez C, Jordá L, et al. Repression of the auxin response pathway increases Arabidopsis susceptibility to necrotrophic fungi. Mol Plant. 2008;1(3):496–509.

    Article  CAS  PubMed  Google Scholar 

  59. Loake G, Grant M. Salicylic acid in plant defence-the players and protagonists. Curr Opin Plant Biol. 2007;10(5):466–72.

    Article  CAS  PubMed  Google Scholar 

  60. Lohse M, Nagel A, Herter T, May P, Schroda M, Zrenner R, Tohge T, et al. Mercator: a fast and simple web server for genome scale functional annotation of plant sequence data. Plant Cell Environ. 2014;37(5):1250–8.

    Article  CAS  PubMed  Google Scholar 

  61. Lorenzo O, Chico JM, Sánchez-Serrano JJ, Solano R. JASMONATE-INSENSITIVE1 encodes a MYC transcription factor essential to discriminate between different jasmonate-regulated defense responses in Arabidopsis. Plant Cell. 2004;16:1938–50.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  62. Lorenzo O, Piqueras R, Sánchez-Serrano JJ, Solano R. ETHYLENE RESPONSE FACTOR1 integrates signals from ethylene and jasmonate pathways in plant defense. Plant Cell. 2003;15:165–78.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  63. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15(12):1–21.

    Article  CAS  Google Scholar 

  64. Malonek S, Rojas MC, Hedden P, Gaskin P, Hopkins P, Tudzynski B. Functional characterization of two cytochrome P450 monooxygenase genes, P450–1 and P450–4, of the gibberellic acid gene cluster in Fusanum proliferatum (Gibberella fujikuroi MP-D). Appl Environ Microbiol. 2005;71(3):1462–72.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  65. Maor R, Puyesky M, Horwitz BA, Sharon A. Use of green fluorescent protein (GFP) for studying development and fungal-plant interaction in Cochliobolus heterostrophus. Mycol Res. 1998;102:491–6.

    Article  Google Scholar 

  66. Meyer FE, Shuey LS, Naidoo S, Mamni T, Berger DK, Myburg AA, van den Berg N, et al. Dual RNA-sequencing of Eucalyptus nitens during Phytophthora cinnamomi challenge reveals pathogen and host factors influencing compatibility. Front Plant Sci. 2016;7:1–15.

    Article  Google Scholar 

  67. Miya A, Albert P, Shinya T, Desaki Y, Ichimura K, Shirasu K, Narusaka Y, et al. CERK1, a LysM receptor kinase, is essential for chitin elicitor signaling in Arabidopsis. Proc Natl Acad Sci U S A. 2007;104(49):1–6.

    Article  Google Scholar 

  68. Montesano M, Brader G, Palva ET. Pathogen derived elicitors: searching for receptors in plants. Mol Plant Pathol. 2003;4(1):73–9.

    Article  CAS  PubMed  Google Scholar 

  69. Mullins E, Chen X, Romaine P, Raina R, Geiser D, Kang S. Agrobacterium-mediated transformation of Fusarium oxysporum: an efficient tool for insertional mutagenesis and gene transfer. Phytopathology. 2001;91:173–80.

    Article  CAS  PubMed  Google Scholar 

  70. Muñoz-Adalia EJ, Fernández M, Wingfield BD, Diez JJ. In silico annotation of five candidate genes associated with pathogenicity in Fusarium circinatum. For Pathol. 2018;48:e12417.

    Article  Google Scholar 

  71. Nagy NE, Fossdal CG, Krokene P, Krekling T, Lönneborg A, Solheim H. Induced responses to pathogen infection in Norway spruce phloem: changes in polyphenolic parenchyma cells, chalcone synthase transcript levels and peroxidase activity. Tree Physiol. 2004;24(5):505–15.

    Article  CAS  PubMed  Google Scholar 

  72. Naidoo S, Visser EA, Zwart L, Du Toit Y, Bhadauria V, Shuey LS. Dual RNA-seq to elucidate the plant– pathogen duel. Curr Issues Mol Biol. 2018;27:127–42.

    Article  PubMed  Google Scholar 

  73. Navarro L, Dunoyer P, Jay F, Arnold B, Dharmasiri N, Estelle M, et al. A plant miRNA contributes to antibacterial resistance by repressing auxin signaling. Science. 2006;312(5772):436–9.

    Article  CAS  PubMed  Google Scholar 

  74. Niño-Sánchez J, Casado-Del Castillo V, Tello V, De Vega-Bartol JJ, Ramos B, Sukno SA, Díaz Mínguez JM. The FTF gene family regulates virulence and expression of SIX effectors in Fusarium oxysporum. Mol Plant Pathol. 2016;17(7):1124–39.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  75. Niño-Sánchez J, Tello V, Casado-del Castillo V, Thon MR, Benito EP, Díaz-Mínguez JM. Gene expression patterns and dynamics of the colonization of common bean (Phaseolus vulgaris L.) by highly virulent and weakly virulent strains of Fusarium oxysporum. Front Microbiol. 2015;6:234.

    PubMed  PubMed Central  Google Scholar 

  76. Norton G, Pappusamy A, Yusof F, Pujade-Renaud V, Perkins M, Griffiths D, Jones H. Characterisation of recombinant Hevea brasiliensis allene oxide synthase: effects of cycloxygenase inhibitors, lipoxygenase inhibitors and salicylates on enzyme activity. Plant Physiol Biochem. 2007;45(2):129–38.

    Article  CAS  PubMed  Google Scholar 

  77. Pieterse CMJ, Van der Does D, Zamioudis C, Leon-Reyes A, Van Wees SCM. Hormonal modulation of plant immunity. Annu Rev Cell Dev Biol. 2012;28(1):489–521.

    Article  CAS  PubMed  Google Scholar 

  78. Qin G, Gu H, Zhao Y, Ma Z, Shi G, Yang Y, Pichersky E, et al. An indole-3-acetic acid carboxyl methyltransferase regulates Arabidopsis leaf development. Plant Cell. 2005;17(10):2693–704.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  79. Rakwal R, Yang G, Komatsu S. Chitinase induced by jasmonic acid, methyl jasmonate, ethylene and protein phosphatase inhibitors in rice. Mol Biol Rep. 2004;31:113–9.

    Article  CAS  PubMed  Google Scholar 

  80. Ramos B, Alves-Santos FM, García-Sánchez MA, Martín-Rodrigues N, Eslava AP, Díaz-Mínguez JM. The gene coding for a new transcription factor (ftf1) of Fusarium oxysporum is only expressed during infection of common bean. Fungal Genet Biol. 2007;44:864–76.

    Article  CAS  PubMed  Google Scholar 

  81. Ren Y-Y, West CA. Elicitation of diterpene biosynthesis in rice (Oryza sativa L.) by chitin. Plant Physiol. 1992;99(3):1169–78.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  82. Richard S, Lapointe G, Rutledge RG, Séguin A. Induction of chalcone synthase expression in white spruce by wounding and jasmonate. Plant Cell Physiol. 2000;41(8):982–7.

    Article  CAS  PubMed  Google Scholar 

  83. Robertson G, Schein J, Chiu R, Corbett R, Field M, Jackman SD, Mungall K, et al. De novo assembly and analysis of RNA-seq data. Nat Methods. 2010;7(11):909–12.

    Article  CAS  PubMed  Google Scholar 

  84. Sanz F, Latour S, Neves M, Bastet E, Pischedda D, Piñeiro G, et al. Industrial applications of Pinus pinaster. Madeira: CIS Madeira, FIBA, AIMMP, CTBA; 2006.

    Google Scholar 

  85. Sasaki Y, Asamizu E, Shibata D, Nakamura Y, Kaneko T, Awai K, Amagai M, et al. Monitoring of methyl jasmonate-responsive genes in Arabidopsis by cDNA macroarray: self-activation of jasmonic acid biosynthesis and crosstalk with other phytohormone signaling pathways. DNA Res. 2001;8(4):153–61.

    Article  CAS  PubMed  Google Scholar 

  86. Schlink K. Identification and characterization of differentially expressed genes from Fagus sylvatica roots after infection with Phytophthora citricola. Plant Cell Rep. 2009;28(5):873–82.

    Article  CAS  PubMed  Google Scholar 

  87. Serra-Varela MJ, Alía R, Pórtoles J, Gonzalo J, Soliño M, Grivet D, Raposo R. Incorporating exposure to pitch canker disease to support management decisions of Pinus pinaster Ait. In the face of cchange. PLoS One. 2017;12(2):1–18.

    Article  CAS  Google Scholar 

  88. Sharma P, Børja D, Stougaard P, Lönneborg A. PR-proteins accumulating in spruce roots infected with a pathogenic Pythium sp. isolate include chitinases, chitonases and β-1,3-glucanases. Physiol Mol Plant Pathol. 1993;43:57–67.

    Article  CAS  Google Scholar 

  89. Simão FA, Waterhouse RM, Ioannidis P, Kriventseva EV, Zdobnov EM. BUSCO: assessing genome assembly and annotation completeness with single-copy orthologs. Bioinformatics. 2015;31(19):3210–2.

    Article  CAS  PubMed  Google Scholar 

  90. Slater GSC, Birney E. Automated generation of heuristics for biological sequence comparison. BMC Bioinformatics. 2005;6:1–11.

    Article  CAS  Google Scholar 

  91. Smith-Unna RD, Boursnell C, Patro R, Hibberd JM, Kelly S, Street D, Brook S, et al. TransRate: reference free quality assessment of de novo transcriptome assemblies. Genome Res. 2016;26(8):21626.

    Article  CAS  Google Scholar 

  92. Soneson C, Love MI, Robinson MD. Differential analyses for RNA-seq: transcript-level estimates improve gene-level inferences. F1000Research. 2016;4(2):1521.

    Article  PubMed Central  Google Scholar 

  93. Surget-Groba Y, Montoya-Burgos JI. Optimization of de novo transcriptome assembly from next-generation sequencing data. Genome Res. 2010.

    Google Scholar 

  94. Tang S, Lomsadze A, Borodovsky M. Identification of protein coding regions in RNA transcripts. Nucleic Acids Res. 2015;43(12):1–10.

    Article  CAS  Google Scholar 

  95. Thaler JS, Humphrey PT, Whiteman NK. Evolution of jasmonate and salicylate signal crosstalk. Trends Plant Sci. 2012;17(5):260–70. Elsevier Ltd.

    Article  CAS  PubMed  Google Scholar 

  96. Thatcher LF, Manners JM, Kazan K. Fusarium oxysporum hijacks COI1-mediated jasmonate signaling to promote disease development in Arabidopsis. Plant J. 2009;58(6):927–39.

    Article  CAS  PubMed  Google Scholar 

  97. Thimm O, Bläsing O, Gibon Y, Nagel A, Meyer S, Krüger P, Selbig J, et al. MAPMAN: a user-driven tool to display genomics data sets onto diagrams of metabolic pathways and other biological processes. Plant J. 2004;37(6):914–39.

    Article  CAS  PubMed  Google Scholar 

  98. Thines B, Katsir L, Melotto M, Niu Y, Mandaokar A, Liu G, Nomura K, et al. JAZ repressor proteins are targets of the SCFCOI1complex during jasmonate signalling. Nature. 2007;448(7154):661–5.

    Article  CAS  PubMed  Google Scholar 

  99. Tiwari SB. Aux/IAA proteins contain a potent transcriptional repression domain. Plant Cell. 2004;16(2):533–43.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  100. Treutter D. Significance of flavonoids in plant resistance: a review. Environ Chem Lett. 2006;4(3):147–57.

    Article  CAS  Google Scholar 

  101. Troncoso C, González X, Bömke C, Tudzynski B, Gong F, Hedden P, Rojas MC. Gibberellin biosynthesis and gibberellin oxidase activities in Fusarium sacchari, Fusarium konzum and Fusarium subglutinans strains. Phytochemistry. 2010;71(11–12):1322–31. Elsevier Ltd.

    Article  CAS  PubMed  Google Scholar 

  102. Truman W, Bennett MH, Kubigsteltig I, Turnbull C, Grant M. Arabidopsis systemic immunity uses conserved defense signaling pathways and is mediated by jasmonates. Proc Natl Acad Sci. 2007;104(3):1075–80.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  103. Tsavkelova E, Oeser B, Oren-Young L, Israeli M, Sasson Y, Tudzynski B, Sharon A. Identification and functional characterization of indole-3-acetamide-mediated IAA biosynthesis in plant-associated Fusarium species. Fungal Genet Biol. 2012;49(1):48–57. Elsevier Inc.

    Article  CAS  PubMed  Google Scholar 

  104. Urban M, Cuzick A, Rutherford K, Irvine A, Pedro H, Pant R, Sadanadan V, et al. PHI-base: a new interface and further additions for the multi-species pathogen-host interactions database. Nucleic Acids Res. 2017;45(D1):D604–10.

    Article  CAS  PubMed  Google Scholar 

  105. Van Bockhaven J, Spíchal L, Novák O, Strnad M, Asano T, Kikuchi S, Höfte M, et al. Silicon induces resistance to the brown spot fungus Cochliobolus miyabeanus by preventing the pathogen from hijacking the rice ethylene pathway. New Phytol. 2015;206(2):761–73.

    Article  CAS  PubMed  Google Scholar 

  106. Van der Does D, Leon-Reyes A, Koornneef A, Van Verk MC, Rodenburg N, Pauwels L, Goossens A, et al. Salicylic acid suppresses jasmonic acid signaling downstream of SCFCOI1-JAZ by targeting GCC promoter motifs via transcription factor ORA59. Plant Cell. 2013;25(2):744–61.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  107. Van Loon LC, Van Strien EA. The families of pathogenesis-related proteins, their activities, and comparative analysis of PR-1 type proteins. Physiol Mol Plant Pathol. 1999;55:85–97.

    Article  Google Scholar 

  108. Visser EA, Wegrzyn JL, Myburg AA, Naidoo S. Defence transcriptome assembly and pathogenesis related gene family analysis in Pinus tecunumanii (low elevation). BMC Genomics. 2018;19(1):1–13. BMC Genomics.

    Article  CAS  Google Scholar 

  109. Visser EA, Wegrzyn JL, Steenkamp ET, Myburg AA, Naidoo S. Dual RNA-Seq analysis of the pine-Fusarium circinatum interaction in resistant (Pinus tecunumanii) and susceptible (Pinus patula) hosts. Microorganisms. 2019;7(9):315.

    Article  PubMed Central  Google Scholar 

  110. Visser EA, Wegrzyn JL, Steenkmap ET, Myburg AA, Naidoo S. Combined de novo and genome guided assembly and annotation of the Pinus patula juvenile shoot transcriptome. BMC Genomics. 2015;16(1):1057. BMC Genomics.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  111. Vivas M, Martín JA, Gil L, Solla A. Evaluating methyl jasmonate for induction of resistance to Fusarium oxysporum, F circinatum and Ophiostoma novo-ulmi. For Syst. 2012;21(2):289–99.

    Article  Google Scholar 

  112. Wan J, Zhang X-C, Neece D, Ramonell KM, Clough S, Kim S-Y, Stacey MG, et al. A LysM receptor-like kinase plays a critical role in chitin signaling and fungal resistance in Arabidopsis. Plant Cell. 2008;20(2):471–81.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  113. Wang X, Tang C, Deng L, Cai G, Liu X, Liu B, Han Q, et al. Characterization of a pathogenesis-related thaumatin-like protein gene TaPR5 from wheat induced by stripe rust fungus. Physiol Plant. 2010;139(1):27–38.

    Article  CAS  PubMed  Google Scholar 

  114. Wegrzyn JL, Liechty JD, Stevens KA, Wu LS, Loopstra CA, Vasquez-Gross HA, Dougherty WM, et al. Unique features of the loblolly pine (Pinus taeda L.) megagenome revealed through sequence annotation. Genetics. 2014;196(3):891–909.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  115. Wegrzyn JL, Lin BY, Zieve JJ, Dougherty WM, Martínez-García PJ, Koriabine M, Holtz-Morris A, et al. Insights into the loblolly pine genome: characterization of BAC and fosmid sequences. PLoS One. 2013;8(9):e72439.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  116. Westermann AJ, Gorski SA, Vogel J. Dual RNA-seq of pathogen and host. Nat Rev Microbiol. 2012;10(9):618–30. Nature Publishing Group.

    Article  CAS  PubMed  Google Scholar 

  117. Wingfield BD, Steenkamp ET, Santana QC, Coetzee MPA, Bam S, Barnes I, Beukes CW, et al. First fungal genome sequence from Africa: a preliminary analysis. S Afr J Sci. 2012;108(1/2):1–9.

    Article  CAS  Google Scholar 

  118. Wingfield MJ, Hammerbacher A, Ganley RJ, Steenkamp ET, Gordon TR, Wingfield BD, Coutinho TA. Pitch canker caused by Fusarium circinatum - a growing threat to pine plantations and forests worldwide. Australas Plant Pathol. 2008;37(4):319–34.

    Article  Google Scholar 

  119. Yin C, Park J-J, Gang DR, Hulbert SH. Characterization of a tryptophan 2-monooxygenase gene from Puccinia graminis f. sp. tritici involved in auxin biosynthesis and rust pathogenicity. Mol Plant-Microbe Interact. 2014;27(3):227–35.

    Article  CAS  PubMed  Google Scholar 

  120. Young, M. D., Wakefield, M. J. & Smyth, G. K. (2010) Goseq : gene ontology testing for RNA-seq datasets Reading data. Manual 1–21. Retrieved from

    Google Scholar 

  121. Zhang Z, Schwartz S, Wagner L, Miller W. A greedy algorithm for aligning DNA sequences. J Comput Biol. 2000;7(1–2):203–14.

    Article  CAS  PubMed  Google Scholar 

  122. Zhao QY, Wang Y, Kong YM, Luo D, Li X, Hao P. Optimizing de novo transcriptome assembly from short-read RNA-Seq data: a comparative study. BMC Bioinformatics. 2011;12(Suppl 14):S2.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  123. Zhu X, Soliman A, Islam MR, Adam LR, Daayf F. Verticillium dahliae’s isochorismatase hydrolase is a virulence factor that contributes to interference with potato’s salicylate and jasmonate defense signaling. Front Plant Sci. 2017;8:1–15.

    Article  Google Scholar 

Download references


We acknowledge Maite Morales Clemente for her excellent technical assistance. The authors further acknowledge the Centre for High Performance Computing (CHPC), South Africa, for providing computational resources to this research project. We thank J.M Diaz-Minguez and V. Tello for their help transforming F. circinatum-GFP.


Laura Hernández was supported by a fellowship from INIA (FPI-INIA) and additional funding for a Short-Term Scientific Mission in the Department of Biochemistry, Genetics and Microbiology, University of Pretoria, Pretoria, South Africa, was provided by Pinestrength Cost Action (FP1406). Financial support for this research was provided by project RTA 2017–00063-C04–01 (Programa Estatal I + D + i, INIA, Spain). EAV was supported through the Technology Innovation Agency (TIA) South Africa, Forest Molecular Genetics Cluster Program. SN was supported by the National Research Foundation (NRF) of South Africa, Y-rated grant program. Opinions expressed, and conclusion arrived at are those of the author(s) and are not necessarily to be attributed to the NRF. Funding bodies had no involvement with study design, data collection, analysis and interpretation or preparations of this manuscript.

Author information

Authors and Affiliations



LHE participate in conception and design of the study, performed the experimental work, data analysis and drafted the manuscript. EAV assisted with computational data analysis, biological interpretation and revised the manuscript. EI participated in the experimental design. RR and SN participated in the conception and experimental design of the study, assisted in biological interpretation and with critical evaluation of the manuscript. All authors have read and approved the final version of the manuscript.

Corresponding author

Correspondence to Sanushka Naidoo.

Ethics declarations

Ethics approval and consent to participate

Not applicable. Pine seedlings used in this study were purchased from ‘Eskalmendi’ nursery (Alava, Spain). No field permissions were necessary to collect the plant samples.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary information

Additional file 1

Symptoms at the shoot tip of inoculated (left side) and mock-inoculated (right side) Pinus pinaster seedlings by the end of the experiment (33 dpi).

Additional file 2.

Statistics for each TransABySS and Trinity assembly. N seq: number of transcripts; N bases: number of bases; Mean length: mean length of the transcripts; N50: N50 value; Ns: number of unknown bases; % GC: guanine and cytosine content; trinity-N: in silico normalized trinity assembly; trinity-nN: non-normalized trinity assemblies. * Best quality preliminary assemblies selected to generate the final assembly.

Additional file 3.

Comparative statistics between normalized (Norm) and non-normalized (N-norm) Trinity preliminary assemblies. Kmer value; % of mapped fragments; % of good mapping; AS: assembly score; OP: optimal score; OC: optimal cutoff; Number of good contigs; % good contigs.

Additional file 4.

BUSCO analysis against the embryophyta lineage database comparing the last Pinus de novo transcriptomes published. P. patula v1.0 [110]; P. patula v2.0 and P. tecunumanii [108].

Additional file 5.

Pinus pinaster de novo transcriptome annotation.

Additional file 6.

Pinus pinaster de novo transcriptome annotation by Mercator tool.

Additional file 7.

mapped reads for each species. Number of differential expressed (DE) genes for Pinus pinaster and DE genes for Fusarium circinatum at each time point in inoculated samples (FDR < 0.05; |log2(Fold Change)| > 0.5). Ppin: P. pinaster; Fcir: F. circinatum;HC: high confident.

Additional file 8.

Principal component analyses (PCA) for Pinus pinaster (above) and Fusarium circinatum (below) rlog data of the differential expression gene analysis (DESeq2). In red: mock-inoculated samples; in blue: inoculated samples at 3 dpi; in green: inoculated samples at 5 dpi; in yellow: inoculated samples at 10 dpi.

Additional file 9.

Clustering of Pinus pinaster and Fusarium circinatum differential expressed (DE) genes. For each cluster with gene ontology (GO) enriched terms, number of genes and percentage for genes are indicated.

Additional file 10.

Significantly enriched GO terms identified from Pinus pinaster genes in each cluster.

Additional file 11:

Phytohormone related differentially expressed (DE) genes in Pinus pinaster.

Additional file 12:

Pathogenesis related (PR) genes differentially expressed (DE) in Pinus pinaster.

Additional file 13:

Significantly enriched GO terms identified from high confidence expressed Fusarium circinatum genes.

Additional file 14:

Hormone related differential expressed (DE) genes in Fusarium circinatum.

Additional file 15:

Fusarium circinatum DE genes related to hormone production with hits in the Pathogen Host Interaction (PHI) database.

Additional file 16:.

RNA-seq data statistics for each sample at each time point, before and after filtering and trimming. Dpi: days post-inoculation; BR: biological replicate, RIN: RNA Integrity Number; Q 30: Phred quality score 30.

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (, which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Hernandez-Escribano, L., Visser, E.A., Iturritxa, E. et al. The transcriptome of Pinus pinaster under Fusarium circinatum challenge. BMC Genomics 21, 28 (2020).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: