Sequencing data and diversity estimates
From a total of 30 pine samples, only 27 were considered for further analysis due to low number of reads. From these, a total of 1,605,754 high-quality reads were obtained in the three sampling sites, after denoising and quality filtering. From the PWN dataset, a total of 13 PWN samples were analyzed, resulting in 778,959 high-quality reads. In both datasets, the sequencing depth was sufficient to describe the diversity and richness of fungal communities as indicated by the rarefaction curves’ plateau (Supplemental Figure S1a and S1b). The estimated diversity indexes for both datasets are presented in Tables 1 and 2.
Overall, samples from both datasets presented a Good’s coverage index (0.99-1.00), corroborating the results from the rarefaction curve analysis (Supplemental Figure S1). For all collection sites in the P. pinaster dataset, samples from “PWN-infected” trees presented a lower number of ASVs (C. Lezírias, 22–42 ASVs; Seia, 18–45 ASVs; and Tróia, 59–162 ASVs), in contrast with “non-infected” trees that were richer with more ASVs per sample (C. Lezírias, 124–249 ASV; Seia, 74–206 ASV; and Tróia, 80–256) (Table 1 and Supplemental Table S1). Likewise, fungal communities from “PWN-infected” trees presented lower diversity (0.70–0.83) than “non-infected” ones, as indicated by the Simpson’s and Shannon’s diversity indexes. For “PWN-infected” condition, statistical differences (p < 0.05) were denoted between locations (C. Lezírias, Seia, and Tróia) in almost all indexes, with exception of Simpson’s (Table 1). Quite the opposite, for the “non-infected” condition, only Shannon and Pielou’s indexes were significant different (p < 0.05) between C. Lezírias and the other two locations (Table 1). The richness of the fungal communities associated with B. xylophilus’ (Table 2), assessed by the number of ASVs was statistically higher (p < 0.05) in C. Lezírias (140 ± 18.3) than in Seia (72.5 ± 6.5) and Tróia (84.5 ± 21.9). However, in terms of diversity, the results followed an opposite trend, with fungal communities from C. Lezírias and Seia being more similar to each other than those from Tróia.
Following our previous results based on cultural-dependent isolation of fungal communities24, we hypothesized that “PWN-infected” and “non-infected” mycobiome differ not only due to the presence of the PWN and disease progression but also as a result of the variations in biogeographic patterns aligned with the timing of disease introduction in the different collection sites. To test these hypothesis, beta-diversity was statistically analyzed (PERMANOVA) (Table 3).
Our premise that fungal communities difference correlate with the presence of PWN (p = 0.001), as well as with the collection site (p = 0.001) were corroborated (Table 3). For the PWN-dataset, statistical differences were observed between collection sites (PERMANOVA, Pseudo-F = 4.199, p = 0.001).
Mycobiome composition, abundance and structure of Pinus pinaster
The composition and abundance of the fungal communities of Pinus pinaster is presented in Fig. 1. Within Fungi kingdom, sequences were assigned onto 1067 ASVs within Fungi kingdom, compromising 5 phyla, 77 orders and 218 genera. Overall, only 141 ASVs were shared between “PWN-infected” and “non-infected” trees, regardless the sampling location (Fig. 1a). Supporting the results obtained from the diversity and richness estimators, the “non-infected” trees presented unique 775 ASVs, while “PWN-infected” trees presented only 151 ASVs (Fig. 1a). The fungal composition of P. pinaster differed between “PWN-infected” and “non-infected” trees, with communities from C. Lezírias and Seia being more similar to each other than to those from Tróia (Fig. 1b). The most abundant phyla, in both “non-infected” and “PWN-infected” trees, were Ascomycota (56.2–82.4%) and Basidiomycota (16.4–43.2%). In “PWN-infected” trees, for both Seia and C. Lezírias, the most common orders for were Saccharomycetales (Ascomycota, 54.4–62.7%), followed by Ophiostomatales (Ascomycota, 12.9–16.0%) and Atratiellales (Basidiomycota, 13–14,4%) (Fig. 1b). In Tróia, the most predominant orders were Hymenochaetales (Basidiomycota, 21.8%), Coniochaetales (Ascomycota, 15.0%), Saccharomycetales (Ascomycota, 11.4%) and Ophiostomatales (Ascomycota, 3.63%). The most abundant genera, in Saccharomycetales order, was Nakazawaea (1.46% in Tróia; 29.5% in Seia; and 30% in C. Lezírias), while in Ophiostomatales order, the most common genera were Ophiostoma (0.83% in Tróia; 5.65% in Seia; and 2.91% in C. Lezírias), Leptographium (0.28% in Tróia; 2.37% in Seia; and 6.37% in C. Lezírias); and Ceratocystiopsis (0.90% in Tróia; 3.12% in Seia; and 4.59% in C. Lezírias). In “non-infected” trees, the most abundant orders were for Tróia, unclassified phylum Ascomycota (un_p_Ascomycota; 24.5%), Capnodiales (14.3%) and Hymenochaetales (11.9%); for Seia, Pleosporales (38.5%), Russulales (14.9%) and Helotiales (12.1%); and for C. Lezírias, unclassified phylum Ascomycota (29.4%), Capnodiales (21.9%), and Tremellales (8%). The orders Saccharomycetales and Ophiostomatales were not detected or less than 5% prevalent.
Composition and abundance of fungal communities of Pinus pinaster: (a) Venn diagram representing the distribution of ASVs across conditions (“non-infected” and “PWN-infected” trees); and (b) relative abundance of the top 20 orders for each condition (“non-infected” and “PWN-infected” trees) and location (C. Lezírias, Seia and Tróia).
A principal co-ordinates analysis (PCoA) was performed to ascertain fungal communities’ structure in terms of sampling location and the presence of the PWN (Fig. 2). The first axis (PCoA 1) explains 21.33% of the variation between samples and clearly separates fungal communities of “non-infected” trees from the “PWN-infected” trees. The second axis (PCoA2) explains 12.6% of the variation in the distribution. In fact, among “non-infected” trees, fungal communities tend to cluster by sampling location, being ASV0021 and ASV0293 (both Pleosporales order, Didymellaceae family, unidentified genus) related with the distribution of Seia while ASV0626 (Capnodiales order) related with Tróia and C. Lezírias. Among “PWN-infected” trees, all samples tend to group regardless the sampling location, being the ASV0014 (genus Proceropycnis) and ASV0131 (genus Nakazawaea) related with their distribution. Together, these axes explain 33.93% of the variability observed between conditions and locations.
Principal coordinates analysis (PCoA) of fungal communities of Pinus pinaster using Bray-Curtis distance matrix. The followings ASVs are related with communities’ distribution: ASV0014 (Basidiomycota, Atractiellomycetes, Atractiellales, Hoehnelomycetaceae, Proceropycnis); ASV0021 (Ascomycota, Dothideomycetes, Pleosporales, Didymellaceae); ASV0131 (Ascomycota, Saccharomycetes, Saccharomycetales, Pichiaceae, Nakazawaea); ASV0293 (Ascomycota, Dothideomycetes, Pleosporales, Didymellaceae); ASV0626 (Ascomycota, Dothideomycetes, Capnodiales).
To validate previous results, a SIMPER analysis was performed to identify the ASVs that represented each condition and location (Supplemental Table S2). Considering only these two conditions, the average similarity within PWN-infected trees was 20.06%, while within “non-infected” trees was 18.96%. These results reflect the heterogeneity of the locations and the trees sampled. More importantly, the average dissimilarity between-condition was 94.59% that emphasizes the differences between PWN-infected and non-infected communities. The top ASVs contributors for this dissimilarity were ASV0014, ASV0131, also assigned in the PCoA analysis, and ASV0502 (Ascomycota, Saccharomycetes, Saccharomycetales, Pichiaceae, Nakazawaea, unidentified species) (Supplemental Table S2).
Occurrence of PWN in the sampled trees
The detection of PWN for each tree sampled was determined by 48 h-extraction from wood pellets, followed by nematode counting24. Intriguing, using ITS2 metagenomics, the fungal community of sample T9 from “non-infected” trees of Tróia showed a composition/abundance closer to “PWN-infected” trees and grouped within this condition, according with the PCoA analysis. To validate PWN classes previously considered, total DNA from each pine sample sequenced was also used for qPCR with a specific probe for B. xylophilus detection (Table 4). Overall, the qPCR results were in accordance with previous detection. For most of the samples, the higher PWN classes (e.g., higher number of PWN extracted) resulted on lower CT values, indicating early detection in qPCR amplification curve. However, the CT value for sample T9 (19.44) was unexpected since PWN class was 0 (e.g., no nematodes extracted), suggesting a positive detection of PWN. Thus, clustering of T9 within the “PWN-infected” trees was supported.
Mycobiome communities’ composition, abundance and structure of Bursaphelenchus xylophilus
From PWN dataset, all sequences were assigned onto 616 ASVs and classified in 5 phyla, 57 orders, and 147 genera. From these, only 58 ASVs were shared between all three locations (Fig. 3a). The main phyla describing the fungal communities of B. xylophilus are Ascomycota (56.18–95.69%) and Basidiomycota (2.31–39.81%). The most predominant orders were different according to sampling location: in Tróia, Orbiliales (Ascomycota, 51.1%), Pezizomycotina_ord_Incertae_sedis (Ascomycota, 27.97%) and Saccharomycetales (Ascomycota, 12.01%); in Seia, Sporidiobolales (Basidiomycota, 34.38%), Pleosporales (Ascomycota, 25.1%), Saccharomycetales (Ascomycota, 15,34%); and, in C. Lezírias, Saccharomycetales (Ascomycota, 24.11%) Spizellomycetales (Chytridiomycota, 12.3%) and Capnodiales (Ascomycota, 9.31%) (Fig. 3a). The order Ophiostomatales was also detected in all samples, but its abundance was less than 1%. On genus level, the most abundant and identified genera, according to sampling location, were: in Tróia, Ciliophora (Ascomycota, 28.0%) and Arthrobotrys (Ascomycota, 15.6%); in Seia, Rhodosporidiobolus (Basidiomycota, 30.9%), Alternaria (Ascomycota, 6.6%) and Groenewaldozyma (Ascomycota, 5.7%); and, in C. Lezírias, Cladosporium (Ascomycota, 7.7%), Nakazawaea (Ascomycota, 6.6%) and Fusarium (Ascomycota, 4.6%).
Composition and abundance of fungal communities of Bursaphelenchus xylophilus: (a) Venn diagram representing the distribution of ASVs across conditions; and (b) relative abundance of the top 20 orders for each location (Lezírias, Seia and Tróia).
The PCoA analysis for the B. xylophilus fungal communities is presented in Fig. 4. The first axis (PCoA1) explains 32.29% of the variation and clearly separates Tróia communities from Seia and C. Lezírias, while the second axis, with 14.92% of the variation, separates communities from C. Lezírias and Seia. The spatial distribution of Tróia’s fungal communities was driven by the ASV0140 (Orbiliales order Orbiliaceae family), the ASV0066 (genus Ciliophora) and ASV0458 (Arthrobotrys cladodes). Fungal communities from Seia were influenced by the abundance of ASV0543 (Pseudopithomyces genus) and ASV0606 (Pleosporales order, Didymellaceae family). Together both axes explain 47.21% of the variability observed between locations.
Principal coordinates analysis (PCoA) of fungal communities of Bursaphelenchus xylophilus using Bray-Curtis distance matrix. The followings ASVs are related with communities distribution: ASV0066 (Ascomycota, Pezizomycotina_cls_Incertae_sedis, Pezizomycotina_ord_Incertae_sedis, Pezizomycotina_fam_Incertae_sedis, Ciliophora); ASV0140 (Ascomycota, Orbiliomycetes, Orbiliales, Orbiliaceae); ASV0458 (Ascomycota, Orbiliomycetes, Orbiliales, Orbiliaceae, Arthrobotrys, Arthrobotrys cladodes); ASV0543 (Ascomycota, Dothideomycetes, Pleosporales, Didymosphaeriaceae, Pseudopithomyces); and ASV0606 (Ascomycota, Dothideomycetes, Pleosporales, Didymellaceae).
As seen for P. pinaster dataset, the SIMPER analysis is in accordance with the PCoA analysis for the B. xylophilus dataset (Supplemental Table S3). The average similarity within-locations ranged between 29.96% (Tróia) and 43.51% (Seia). The average dissimilarity between-locations were considerably high (71.08–89.46%), highlighting their distribution in the Fig. 4. For instances, ASV0543 and ASV0606 were among the most representative to separate fungal communities of C. Lezírias and Seia (as in the PCoA). In the case of the B. xylophilus’ fungal communities of C. Lezírias and Tróia, the ASV0140, ASV0066 and ASV0458 were the highest contributors, which are also assigned in Fig. 4.
