In May 2019, Scots pine seedlings (ca 1-year-old, 1,447 seedling plants in total) consisting of 85 full-sib (FS) families, 35 half-sib population mix (HSpm) families, and 116 half-sib (HS) families (Supplementary Table 1) were provided by the Natural Resources Institute Finland (Luke). Each family had six plants, three treated as control (wounds without fungal infection) and three as infected (wounds plus fungal infection). Seedlings were potted using FPM 420 peat substrate (Kekkilä Professional, Finland), and those from different treatments or family types were placed on separate greenhouse tables. Seedlings within the same treatment and family type were randomly arranged. HSpm was originally composed of 50 families from Poland, southern Sweden (Gullabo), southern Finland (Uusikaupunki), central Finland (Kälviä), and northern Finland (Kolari), with 10 families per region. In total, 35 families survived the transfer to the greenhouse, 2, 7, 7, 9, and 10 families from Poland, southern Sweden, southern, central and northern Finland, respectively.

After a 2-month acclimatization in the greenhouse, needle tissues of each plant were sampled before inoculation and stored at −80°C for further use. One H. annosum isolate (homokaryon, 03007) previously shown to induce strong necrotic lesions was cultured at 20°C on malt extract–sawdust agar plates (MEA-S) containing 2% Scots pine sawdust, 2% malt extract, and 1.8% agar for 2 weeks before inoculation. In July 2019, inoculations were conducted on the stem at ca 15 cm above the greenhouse table according to the following steps: surface sterilization of the stem with 70% ethanol, drilling a hole with the puncher of 2 mm in diameter, and filling the punched hole with MEA-S agar plugs pre-colonized by H. annosum. Control treatments received mock inoculation after wounding. Both treatments were sealed with Parafilm M (Heathrow Scientific, USA). At 3 months post-inoculation, the height of each plant was measured from the bottom to the top branches. At 4 months post-inoculation, sapwood lesions of all plants were measured. Lesion size was the sum of vertical and horizontal lesions as described by Keriö et al. (2015). Lesion area was measured using Fiji-ImageJ 1.51 software (Rueden et al., 2017). Lesion area (family mean) of the control and infected plants in each family type were compared by a two-tailed unpaired t-test. Lesion area of HSpm infected plants from three regions in Finland was compared using one-way analysis of variance (ANOVA) followed by Tukey's multiple comparisons at 0.05 confidence level. Pearson correlation coefficients and p values between lesion size, lesion area, and height based on family means by treatments were computed using the corr. test function in the package psych in R 4.1.0 (R Core Team, 2021; Revelle, 2021).

Heritability of lesion area was estimated at an individual (hHS2, hFS2) and family-mean basis (hHSfam2, hFSfam2) separately in the HS and FS data sets using the following equations:

where σGCA2 is the variance of the general combining ability (GCA), σSCA2 is the variance due to the specific combining ability (SCA), σE2 is the error (residual) variance, and n_H is the harmonic mean of the ramets per family (in this case, n_H = 6.11 and n_H = 5.99 for HS and FS, respectively). The variance components were derived from two mixed models where Treatment was a fixed factor and female (HS) or female and male parents (FS) were the random factors. The GCA referred either to the effect of the female parent (HS) or to the average effect of female and male parents (FS). The SCA was an interaction of a specific female–male combination (only estimable from the full-sib data). The analysis was performed using the ASReml software (Gilmour et al., 2015).

Terpene compounds of needles sampled before inoculation from 12 resistant and 12 susceptible families were quantified as described earlier (Liu et al., 2022). Approximately 100 mg fresh weight of needle tissues was ground in liquid nitrogen and was then extracted in 2 mL of n-hexane at room temperature for 2 h and washed twice with 2 mL n-hexane. The compound 1-chloro-octane (70.1 μg) was used as an internal standard. The extracts were analyzed using a 6890A gas chromatograph (Agilent Technologies, USA) equipped with a HP-5MS UI column (30 m × 0.25 mm × 0.25 μm, Agilent J&W Scientific, USA) and connected to a 5,973 inert mass selective detector (Agilent Technologies, USA). Helium was used as the carrier gas, and linear velocity was about 40 cm/s. The pulsed splitless sampling technique was used with the pulse pressure 25 psi and duration of 0.5 min, and 1 μl was injected. The initial oven setpoint temperature of 50°C was held for 1 min and the column temperature was increased from 50 to 115°C at 5°C/min, then to 280°C at 15°C/min and held for 10 min. Mass numbers from m/z 33 to 350 were recorded. Authentic standards, ChemStation software, as well as Wiley275 and Nist17 libraries were used to identify and quantify monoterpenes and sesquiterpenes. Compounds for which authentic standards were not available were quantified using a compound with a similar chemical structure. For each terpene compound, concentrations μg/mg fresh weight of needle tissues were shown as the mean of 36 resistant or 36 susceptible seedlings (12 resistant and 12 susceptible families, 3 seedlings per family), and the comparison between the two groups was assessed by two-tailed unpaired t-test. Pearson correlation coefficients with statistical significances between lesion area and terpene concentrations were computed using the same method as described above.

Total RNA was extracted from ca 100 mg of ground needle tissues sampled prior to the inoculation using the method previously described by Chang et al. (1993). The ground sample was transferred to a sterile 2 ml Eppendorf tube, and 900 μl extraction buffer (preheated at 65°C) and 9 μl DTT (1 mol/L) were added. The extraction buffer contained 2% CTAB, 4% PVP (K30), 100 mM Tris-HCl, 25 mM EDTA, and 2 M NaCl. The mixture was vortexed and incubated at 65°C for 15 min. Equal volumes of chloroform: isoamyl alcohol (24:1) were added, followed by mixing well, and centrifugation at 10,000 g at room temperature for 10 min. The upper phase was transferred to a new 2 ml Eppendorf tube and then equal volumes of chloroform: isoamyl alcohol (24:1) were added and the above step was repeated. To precipitate RNA, 1/4 volume of 10M LiCl (42.4 g/mol) was added, mixed well, and kept at 4°C overnight. The sample was centrifuged at 10,000 g at 4°C for 30 min. The supernatant was pipetted out and the pellet was washed by adding 100 μl cold 70% ethanol. The centrifugation was repeated for 5 min. The supernatant was pipetted out and the pellet was dried in the sterile hood. The pellet was resuspended in 10–20 μl nuclease-free water. Extracted RNA was quantified by NanoDrop 2000c (Thermo Fisher Scientific, USA) and quality checked with Agilent 2100 bioanalyzer (Agilent Technologies, Germany). Equal amounts (1 μg) of total RNA from three seedlings per family were pooled together as one sample. Samples of 12 resistant and 12 susceptible families (36 resistant and 36 susceptible seedlings) for terpene detection were also used for RNA-seq. RNA-seq was performed with paired-end (150 bp) sequencing using an Illumina NovaSeq 6000 platform at Novogene (UK). RNA-seq reads for 24 samples are available in the NCBI GEO database with the accession number GSE200311. SortMeRNA v4.2.0 was used to remove ribosomal RNAs from raw reads (Kopylova et al., 2012), followed by the trimming of low-quality reads at Q < 20 in 5-base sliding windows and adaptor sequences by Trimmomatic v0.39 with a minimum length of 50 bp (Bolger et al., 2014). Trimmed reads were qualified by FastQC v0.11.8 and MultiQC v1.8. De novo transcriptome assembly was constructed with trinity 2.11.0 (Grabherr et al., 2011), with setting normalize max read cov 50 and min kmer cov 2. A trinity assembly with a size of 288 Mb was produced with GC content 42.08% and E90N50 value 1,824 bp. Sample-specific alignment was conducted with bowtie2 (Langmead and Salzberg, 2012). Alignment-based quantification was performed with RNA-Seq expression estimation by Expectation-Maximization (RSEM) (Li and Dewey, 2011).

Differential expression analysis for each family type was carried out using the package DESeq2 in R (Love et al., 2014). Transcripts that were expressed over one count in at least four libraries were kept. Transformed data by variance stabilizing transformation (VST) were input to perform principal component analysis (PCA). To produce a list of differentially expressed transcripts (DETs), cut-offs were set as Benjamini-Hochberg adjusted p value (FDR) < 0.05 and fold change (FC) > 2. Gene ontology (GO) enrichment of DETs was analyzed using the R package clusterProfiler (Wu et al., 2021). GO terms of all Trinity assembly transcripts were annotated through Blast2GO (Conesa et al., 2005). The description file of all GO terms (go-basic.obo) was downloaded from the Gene Ontology database. Significantly enriched terms (FDR < 0.05) were eliminated due to gene overlaps using the function reduce_overlap in the package GOplot in R (Walter et al., 2015). The R package WGCNA was used for co-expression network analysis of HS and HSpm samples (Langfelder and Horvath, 2008). Lowly expressed transcripts with count <1 in all libraries and count <100 in at least eight libraries were removed followed by VST transformation. Gene modules were detected using 1-step network construction function blockwiseModules with proper soft-thresholding power, minimum module size of 100, and merging threshold function of 0.25. Hub genes in the key module were identified using cut-offs of eigengene-based connectivity (KME) ≥ 0.8 and geneTraitSignificance ≥ 0.2.

Necrotic lesions observed in the stem sapwood caused by H. annosum infection were larger than those of the control samples (Figure 1A and Supplementary Table 2). The average lesion area of control and infected seedlings was 3.2 mm² (SE 0.07 mm²) and 6.3 mm² (SE 0.09 mm²), respectively, with significant differences at p < 0.001 for the comparison between the two treatments in each family type. Variations in lesion size, lesion area, and plant height were detected among families and data of the above traits were normally distributed (Figure 2). Lesion area of infected samples ranged from 3.8 to 10.8, 4.0 to 9.3, and 2.3 to 7.9 mm² for FS, HS, and HSpm, respectively (Figure 1C and Table 1). Lesion size was strongly correlated with lesion area. Plant height showed a weak to moderate relationship with lesion area, of which coefficients were 0.18 and 0.37 for control and infected samples, respectively (Figure 2). Individual heritability estimates of lesion area of FS and HS were similar, around 0.35, whereas FS represented a larger value on the family-mean basis at 0.53 than HS at 0.38 (Table 2).

Seedlings were considered to be resistant if they had lesion areas of 2.3–4.4 mm² while susceptible seedlings had a lesion area of 6.6–9.9 mm². Based on the family-mean lesion area, 12 resistant and 12 susceptible families of FS, HS, and HSpm were selected for further terpene and RNA-seq analyses (Table 1). Lesions displayed obvious differences between resistant and susceptible seedlings (Figure 1D). Compared with FS and HS, HSpm seedlings showed a smaller lesion area. Families 302, 296, and 300 in HSpm from southern Finland showed enhanced resistance, whereas 307, 303, 311, and 308 from the north were mostly susceptible. Means of lesion area of infected samples increased from the south (3.9 mm²) to central (4.9 mm²) and north (6.3 mm²) parts of Finland (Figure 1B and Supplementary Table 3).

From needle tissues of 36 resistant and 36 susceptible seedlings which were sampled prior to inoculation, 17 monoterpene and 12 sesquiterpene compounds were quantified (Supplementary Table 4). Monoterpenes were the predominant terpene compounds, of which total concentration was ca 20-fold higher than sesquiterpenes. α-Pinene and β-caryophyllene showed the highest concentrations among monoterpenes and sesquiterpenes, respectively. A few monoterpenes were significantly different between resistant and susceptible seedlings (Table 3). 3-Carene was highly abundant in resistant genotypes compared to susceptible ones. By contrast, the susceptible group had highly concentrated tricyclene, β-pinene, camphene, limonene+β-phellandrene, p-cymene, and α-pinene. Contents of these compounds were positively highly correlated with each other but negatively correlated with 3-carene (Supplementary Table 5). Correspondingly, lesion area was negatively correlated with 3-carene content with a coefficient at −0.34, but positively with tricyclene and β-pinene at 0.28 and 0.26, respectively.

Each library obtained 37–54 million paired-end reads, which were de novo-assembled through the Trinity platform with alignment rates of 91.4–98.9%. Although we mapped the reads against an alternative reference genome of Pinus taeda (Pyhäjärvi et al., 2020), the mapping rates were very low for each library. Samples differing in resistance to H. annosum were separated along the y-axis, x-axis, and the diagonal for FS, HS, and HSpm, respectively, in PCA plots (Figure 3).

Differential expression analysis of each family type identified transcripts that were highly expressed either in resistant or susceptible seedlings (Supplementary Table 6). Functional analysis of subsets of DETs showed two significantly enriched GO terms in HSpm and 23 terms in HS (Figure 4). Resistant samples of HSpm were enriched in DETs involved in the flavonol biosynthetic process and response to other organisms. While HS-resistant samples were enriched in auxin metabolic process and linoleate 9S-lipoxygenase activity, other terms such as hormone activity, cellular response to heat, and wax biosynthetic process were mostly prevalent in susceptible genotypes.

Among DETs of each family type, 36–57% were annotated (Supplementary Tables 6, 7). Signaling was primarily mediated by plant hormones such as jasmonic acid, abscisic acid, and gibberellins (Figure 5). Calcineurin B-like (CBL) interacting-protein kinases (CIPKs), RALF-like proteins (RALFLs), and calcium-dependent protein kinases (CDPKs) in susceptible seedlings could play a role to regulate calcium-dependent signaling. Despite differences in resistant capability against H. annosum, all samples had highly expressed transcripts involved in defense encoding disease resistance proteins, pathogenesis-related proteins, and key enzymes in the phenylpropanoid and terpenoid biosynthesis. Susceptible samples showed increases in DETs such as MYB transcription factors (MYBs) and triphosphate tunnel metalloenzyme 2 (TTM2), as well as abundant DETs associated with abiotic stress responses (Figure 5). Susceptible samples of HS had increased levels of 22 transcripts encoding heat shock proteins (HSPs), but negligible transcript levels were observed in other samples.

In addition to functional characterization of DETs, we applied WGCNA to shed light on co-expression networks between transcripts according to their common expression trends across all half-sib (HS and HSpm) samples. After filtering, 20,618 transcripts were used to construct the network to find key modules of highly correlated genes. Samples from the same family type, having a similar lesion area clustered well (Figure 6A), where red means a large lesion and white a small lesion. Using the soft-thresholding power of 6, we identified 49 modules (Supplementary Figures 1A,B). Module size ranged from 123 to 2,448 transcripts. There were 957 transcripts in the gray module that failed to cluster to any module. The dark gray, saddle brown, and pale turquoise modules were significantly correlated with lesion area with coefficients at −0.53, 0.59, and 0.63, respectively (Supplementary Figure 1C). We considered the pale turquoise as the key module since it was closely clustered with the lesion area (Figure 6B). A total of 196 transcripts in the module were highly correlated with the lesion area (Figure 6C). Among them, 72 transcripts were hub genes with high connectivity, and 58 of them were annotated with 20 transcripts for HSPs (Supplementary Table 8), mostly enhanced in the susceptible genotypes. Hub genes are transcripts with high intramodular connectivity, which tend to be strongly correlated with the lesion trait.