Acquisition of fungi from the environment modifies ambrosia beetle mycobiome during invasion

Ambrosia beetle species

We selected two ambrosia beetle species: the exotic X. germanus and the native X. saxesenii. Xylosandrus germanus is a species native to Asia that was first reported in Europe in the 1950s and since then rapidly spread, becoming one of the dominant ambrosia beetles in European forest ecosystems (Galko et al., 2018). Xylosandrus germanus’ fungal mutualist is Ambrosiella grosmanniae (Mayers et al., 2015). Xyleborinus saxesenii is instead a species of Palaearctic origin and its main fungal mutualist is Raffaelea sulfurea, although other fungi have been found in association with this beetle species (Biedermann et al., 2013).

Sampling locations and procedure

Beetles were collected in 2016 in ten forest stands located in the Northeast of Italy (Fig. S1 and Table S1), across two forest habitats: old-growth forests (n = 5) and restored forests (n = 5). With “old-growth forests”, we refer to the remnants of the old oak–hop-hornbeam forest (Quercus spp. and Ostrya carpinifolia Scop.) that covered the vast majority of Veneto and Friuli Venezia Giulia regions after the last ice age. With “restored forests”, we refer to mixed forests that were planted over the last 30 years to restore forests across agricultural landscapes. Both forest habitats are dominated by oak (Quercus spp.), ash (Fraxinus spp.), maple (Acer spp.), and hop-hornbeam (O. carpinifolia). In addition, both forest habitats are present in relatively small patches embedded in an agriculture-dominated landscape (min = 2.65 ha, max = 165.15 ha for old-growth forests; min = 2.37 ha, max = 37.41 ha for restored forests).

Beetles were trapped using green and purple 12-multi-funnel traps (Synergy Semiochemicals, Burnaby, Canada) baited with ultra-high release rate ethanol pouches (99% purity, release rate of 300–400 mg/day at 20 °C, Contech Enterprises). Although ethanol is attractive for a wide range of wood-borers (Miller, 2006), it is also the most commonly used volatile for trapping ambrosia beetles (Reding et al., 2011). The ethanol pouch was always attached to the sixth funnel and hung outside the trap. Traps were set in the understory at about 1.5 m above the ground and were suspended at least 1m from the tree bole. Trap collecting cups were half-filled with 1:1 solution (v/v) of ethylene glycol:water to kill and preserve captured beetles (known as “wet system”) (Steininger et al., 2015). At each trap check, collecting cups were emptied and the solution was renewed. Traps were set up in mid-May and emptied every three weeks until the beginning of August. At each visit, all insects were collected, put in tubes filled with ethanol, and brought to the laboratory where ambrosia beetles were separated from other trapped insects. Each individual was then morphologically identified to species level and kept in separate vials filled with ethanol until they were processed. Then, we retained for the analysis only individuals of the two ambrosia beetle species that were simultaneously collected during the same trapping period and in the same trap. This allowed an intra-trap comparison to test for cross-contamination between beetle species (see ‘Results’, Table S2). Our sampling procedure did create the possibility of microbial cross-contamination among the different insect specimens simultaneously present in the trap collector cup (Viiri, 1997). However, in our previous work we demonstrated that individuals collected in the same trap do not show evidence of cross-contamination (Malacrinò et al., 2017). In an effort to reduce possible environmental contamination, we also sterilized the external surface of the insect body. First, we put each insect in a vial with ddH₂O in a water bath and sonicated them for 1 min. After sonication, we washed each insect by vortexing once in ethanol (100%), twice in sodium hypochlorite (5%), and twice in ddH₂O for 1 min following each wash step. For each ambrosia beetle species, we processed 15 individuals per sampling site (total of 300 individuals).

DNA extractions, libraries preparation and amplicon sequencing

Single individuals were crushed in an extraction buffer (10 mM Tris, 100 mM NaCl, 10 mM EDTA, 0.5% SDS) using three one mm ∅ stainless steel beads per tube, with the aid of a bead mill homogenizer set at 30 Hz for 5 min (TissueLyzer II, Qiagen, UK). The mixture was treated with proteinase K (5Prime GmbH, Germany) following the producer’s instructions. Total DNA was extracted using the MoBio PowerSoil Kit (Mo Bio Laboratories, Inc., Carlsbad, CA, USA) following the manufacturer’s protocol. DNA concentration and purity were assessed with a Nanodrop 2000 spectrophotometer (Thermo Fisher Scientific Inc., USA).

The fungal community associated with each individual was characterized by amplicon sequencing targeting the ITS2 region using gITS7 and ITS4 primers, as previously indicated by Kostovcik et al. (2015). We selected the ITS2 region (Nilsson et al., 2019) to ensure we captured the diversity of fungi with which beetles come in contact during host searching. We are aware that the ITS2 region can lead to an amplification bias for Microascales and Ophiostomatales (Kostovcik et al., 2015), the two orders including the main mutualists of X. germanus and X. saxesenii. Here, however, we are interested in the entire mycobiota, which has been less frequently described and might explain important aspects of beetle ecology. PCR reactions were performed in a total volume of 25 µl, containing about 50 ng of DNA, 0.5 µM of each primer, 1X KAPA HiFi HotStart ReadyMix (KAPA Biosystems, USA) and nuclease-free water. Amplifications were performed in a Mastercycler Ep Gradient S (Eppendorf, Germany) set at 95 °C for 3 min, 98 °C for 30 s, 56 °C for 30 s and 72 °C for 30 s, repeated 30 times, and ended with 10 min of extension at 72 °C. Each amplification was carried out in technical triplicate, including three non-template controls with nuclease-free water instead of DNA. Nuclease-free water (100 µl) was also processed using the same procedure as the experimental samples (from DNA extraction to sequencing) in order to exclude contamination of reagents and instruments. Amplification success was checked by electrophoresis on 1.5% agarose gel stained with GelRed (Biotium Inc., Fremont, CA, USA). Although we did not observe any amplification bands for the negative/non-template control samples, these were processed and sequenced together with experimental samples. PCR products from the same sample were then pooled together, and cleaned using Agencourt AMPure XP kit (Beckman Coulter, Brea, CA, USA) following the producer’s instructions. A further short-run PCR was performed to integrate Illumina i7 and i5 indexes following the producer’s protocol (Nextera XT, Illumina, San Diego, CA, USA), and amplicons were purified again as explained above. Libraries were then quantified with the Invitrogen Qubit HS dsDNA kit (Invitrogen, Carlsbad, CA, USA), normalized to a concentration of 10 ng/µl using nuclease-free water, pooled together and sequenced with an Illumina MiSeq sequencer, using the MiSeq Reagent Kit v3 300PE chemistry (Illumina, San Diego, CA, USA) following the producer’s protocol.

Data analysis

Demultiplexed paired-end reads were merged using the PEAR 0.9.1 algorithm using default parameters (Zhang et al., 2014). Raw data handling was carried out using QIIME 1.9 (Caporaso et al., 2012), quality filtering reads with default parameters, binning Operational Taxonomic Units (OTUs) using open-reference OTU-picking through UCLUST algorithm (97% similarity), and discarding chimeric sequences discovered with USEARCH 6.1 (Edgar, 2010). All non-fungal OTUs were discarded using ITSx (Bengtsson-Palme et al., 2013). Taxonomy assignment was performed using the BLAST method (default parameters) by querying towards a custom database built using all ITS2 reference sequences deposited at NCBI GenBank (accessed on July 2017). R statistical environment v3.5.1 (R Core Team, 2013) plugged with the packages vegan (Dixon, 2003), phyloseq (McMurdie & Holmes, 2013), picante (Kembel et al., 2010) and DESeq2 (Love, Huber & Anders, 2014) was used for data analysis. First, singletons and samples with fewer than 1,000 counts were removed. Data processing resulted in a dataset of 3,634,647 reads clustered into 19,744 OTUs. Then, comparisons of fungal community composition between ambrosia beetle species and between forest habitats within the same beetle species were performed using PERMANOVA analysis (999 permutations stratified at site level) calculated on a UniFrac distance matrix. Non-metric multidimensional scaling (NMDS) procedure was performed to visualize differences in the structure of fungal communities. The diversity of fungal communities was assessed using Chao1 (total diversity) (Chao, 1984), Faith’s phylogenetic diversity (which considers both total diversity and the phylogenetic relationship between taxa within the community) (Faith, 1992), and 1-Simpson (dominance) (Simpson, 1949) indexes. Comparisons were performed using mixed-effects models (one model for each diversity index) with the lmer function under the lme4 R package (Bates et al., 2015) using ambrosia beetle species and forest habitat as factors, and sampling site as a random variable. The package emmeans was used to infer pairwise contrasts within mixed-effects models (FDR corrected). The use of “sampling site” for stratification in PERMANOVA and as a random variable in the mixed-effects model allowed the control of both non-homogeneity in sample number at each trap (Table S2) and potential spatial effects. The differential presence of OTUs between forest habitats and within the same beetle species was assessed using the package DESeq2, by contrasting the two forest types within each species. The association of each fungal genus to a functional guild was performed by searching against the FUNGuild database (Nguyen et al., 2016) and manually curating the results in case of multiple results from the same query.

Fungal communities associated with X. germanus and X. saxesenii

The reconstruction of the fungal communities showed that the exotic ambrosia beetle X. germanus and the native ambrosia beetle X. saxesenii are associated with different fungi (F_1,211 = 10.5; P < 0.001). The absence of cross-contamination was shown by a multiple comparison procedure following PERMANOVA: differences between ambrosia beetle species were found at all sites (P < 0.01 FDR corrected–Table S2). In case of cross-contamination we would expect an overlap of the fungal communities and, thus, no differences.

The fungal communities of both ambrosia beetle species were dominated by unidentified taxa (83.97% in X. germanus and 73.91% in X. saxesenii–Table S3). Instead, we identified 26 genera that include plant pathogens (4.75% in X. germanus and 5% in X. saxesenii) mainly represented by the genus Cladosporium. The rest of the communities were represented by saprotrophs (2.73% in X. germanus and 4.12% in X. saxesenii), yeasts (5.35% in X. germanus and 15.96% in X. saxesenii) (Table S3) and at low relative abundances insect pathogens, mycorrhizal fungi, endophytes and lichen parasites (Table S3). In addition, in both ambrosia beetle species we found sequences that can likely be assigned to the respective main mutualists: Ambrosiella sp. in X. germanus (2.25%) and Raffaelea sp. in X. saxesenii (0.01%).

Effect of forest habitat on fungal communities

Using a PERMANOVA analysis we found that the fungal community segregated by forest habitat both in the exotic X. germanus (F_1,10 = 21.8; P < 0.001—Fig. 1A) and in the native X. saxesenii (F_1,10 = 1.6; P = 0.004—Fig. 1B), although the effect was much more evident in X. germanus (R² = 0.16) than in X. saxesenii (R² = 0.01). A different pattern between X. germanus and X. saxesenii also emerged when looking at the diversity of the fungal communities. For X. germanus, individuals collected in old-growth forests were associated with a richer and more diverse fungal community than those collected in restored forests (P < 0.001 for both Chao1 and phylogenetic diversity—Figs. 2A and 2B, Table S4), whereas these differences were not observed in the native X. saxesenii (P > 0.05 for both Chao1 and phylogenetic diversity—Figs. 2A and 2B, Table S4). On the contrary, in both X. germanus and X. saxesenii the dominance index (1-Simpson) significantly differed between the two forest habitats (P < 0.001—Fig. 2C, Table S4). In particular, we observed a higher dominance index in restored forests compared to old-growth forests for X. germanus, and a higher dominance index in old-growth forests versus restored forests for X. saxesenii.

Comparing the fungal community associated with individuals collected in the two forest habitats, for X. germanus we found 121 differentially abundant OTUs: 4 of them were more abundant in restored than in old-growth forests (1 Ambrosiella sp., 1 Aspergillus sp., 1 Saccharomyces sp. and 1 unidentified, Fig. 3A), whereas 117 were more abundant in old-growth than in restored forests (102 unidentified OTUs, and 15 genera, Fig. 3A). The same analysis on X. saxesenii resulted in 4 differentially abundant OTUs (1 Aureobasidium sp. and 3 unidentified), all of them more abundant in restored than in old-growth forests (Fig. 3B).