Bacterial community associated to the pine wilt disease insect vectors Monochamus galloprovincialis and Monochamus alternatus

Monochamus beetles are the dispersing vectors of the nematode Bursaphelenchus xylophilus, the causative agent of pine wilt disease (PWD). PWD inflicts significant damages in Eurasian pine forests. Symbiotic microorganisms have a large influence in insect survival. The aim of this study was to characterize the bacterial community associated to PWD vectors in Europe and East Asia using a culture-independent approach. Twenty-three Monochamus galloprovincialis were collected in Portugal (two different locations); twelve Monochamus alternatus were collected in Japan. DNA was extracted from the insects’ tracheas for 16S rDNA analysis through denaturing gradient gel electrophoresis and barcoded pyrosequencing. Enterobacteriales, Pseudomonadales, Vibrionales and Oceanospirilales were present in all samples. Enterobacteriaceae was represented by 52.2% of the total number of reads. Twenty-three OTUs were present in all locations. Significant differences existed between the microbiomes of the two insect species while for M. galloprovincialis there were no significant differences between samples from different Portuguese locations. This study presents a detailed description of the bacterial community colonizing the Monochamus insects’ tracheas. Several of the identified bacterial groups were described previously in association with pine trees and B. xylophilus, and their previously described functions suggest that they may play a relevant role in PWD.

Bacterial community composition. Seventeen samples were selected for pyrosequencing, 11 samples of M. galloprovincialis and 6 of M. alternatus (Table 2). Sample selection was based on DGGE analysis by picking several samples from the identified clusters.
A total of 185 275 raw reads were obtained that were filtered. In the dereplication step, 38 136 unique reads were identified. Singletons and chimeras were removed resulting in a final number of 228 OTUs. Twelve OTUs corresponded to eukaryotic DNA and were discarded from the analysis, resulting in a total of 216 OTUs and 159 280 high-quality reads ( not be assigned (Supplementary Table S1). Two of the 210 assigned OTUs were only identified up to phylum level and 11 were only identified up to class level. From the 197 OTUs identified up to order level, 179 OTUs were identified at family level, 86 OTUs were identified at genus level and 16 OTUs were identified at species level (Supplementary Table S1). The rarefaction curve for each insect sample tended to saturation ( Figure S1), indicating that the OTUs detected were representative of the bacterial communities.
Monochamus galloprovincialis bacterial community. A total number of 101 091 reads were obtained for M. galloprovincialis samples corresponding to 57 OTUs assigned to 5 different phyla: Proteobacteria, Firmicutes, Bacteroidetes, Acidobacteria and Actinobacteria. Proteobacteria was the dominant phylum, representing on average 98.4% of the reads in each sample with a standard deviation of ±2.0% (Supplementary  Table S2). Firmicutes were present in 5 M. galloprovincialis samples and Bacteroidetes in 9 samples (out of 11). Firmicutes accounted on average for 0.7 ± 1.2% of the total reads in each sample and Bacteroidetes abundance averaged 0.7 ± 0.8%. Acidobacteria and Actinobacteria were each detected in one sample from Comporta, both with less than 1% of the total number of reads.
Twenty one different bacterial orders were detected in M. galloprovincialis samples (Supplementary Table S3). The most represented orders were Vibrionales (24.9 ± 13.9%), Enterobacteriales (34.9 ± 29.3%), Pseudomonadales (18.2 ± 7.7%) and Oceanospirillales (11.6 ± 7.5%) representing altogether 91.5% of the total reads for the Portuguese host and detected in all samples. However, the relative abundance of each order differs among samples (Fig. 2a). For instance, from a total of 11 samples, Vibrionales was the most abundant order in 6, representing 36.2 ± 3.4% of the total reads in these 6 samples; Enterobacteriales was the most abundant order in 4 samples accounting for 72.1 ± 10.5%; Pseudomonadales was the most abundant order in one sample (MG.PC42) with 34.5% of the total reads (Supplementary Table S3).
The most abundant OTUs in our M. galloprovincialis dataset affiliated with the genus Vibrio (OTU 5, 17.1% of the total reads) and with the family Enterobacteriaceae (OTU 2, 10.7% of the total reads). From these, only OTU 5 (Vibrio) was present in all M. galloprovincialis samples. OTUs that were not only abundant (between 3 and 10% of the total reads) but also present in all M. galloprovincialis samples, affiliated with Pseudomonas (OTU 3), Halomonas (OTU 4), Oceanicola (OTU 9), Vibrio (OTU 28) and with Enterobacteriaceae family (OTU 6) (Fig. 3). PERMANOVA analysis of the OTU abundances confirmed that there were no significant differences between the bacterial communities of insects of the two sampling sites (Comporta and Mortágua) (p > 0.05).
Despite the resemblances of the bacterial communities of both insect species, through NMDS analysis ( Fig. 5) it was possible to observe a separation between M. galloprovincialis samples and the ones from   . Heatmap for the OTUs that represent more than 10% of the total number of reads. The colour gradient goes from light yellow for the lower numbers of reads, to dark green for the higher numbers of reads. The numbers in the third top line represent the last two digits in the sample name (see Table 2).
The number of OTUs was much higher in M. alternatus (199) than in M. galloprovincialis (57) with significant differences in species richness between insects species (P < 0.05 Mann-Whitney test). There were no significant differences among the diversity and equitability indexes between the tracheal bacterial communities from M. alternatus and M. galloprovincialis (p > 0.05, Student's t-test).

Discussion
The goal of this study was to characterize the bacterial community colonizing the trachea of M. galloprovincialis and M. alternatus, two insect vectors of B. xylophilus, the causative agent of PWD 5 . The trachea of the insect vector is the preferential organ for nematode lodging 12 .
The role of bacteria in PWD is a controversial subject 23 but nowadays the microbiome is acknowledged as very important to host fitness and adaptation. In insect-plant interaction, microbes are described as crucial mediators, especially for herbivorous insects 24 . Thus, characterizing the microbiome of the PWD insect-vectors is an important step to better understand the disease mechanism and the role of bacteria in the pathogenic process. Considering the nematode life cycle and its relation with the insect, it is reasonable to place the focus on the microbiome of the trachea 9 .
To deeply characterize the bacterial community of the insects' tracheas, a culture-independent approach was used, combining DGGE analysis and barcoded pyrosequencing. DGGE analysis is described as a costand time-effective technique, suitable to evaluate compositional variation of complex bacterial communities 25 . This technique has been used in combination with massive parallel sequencing in many studies, allowing a global view of the community structure of a high number of samples and the selection of a smaller set for deep characterization 26 .
The results from this study show that the bacterial community of both Monochamus species is dominated by Proteobacteria, mostly represented by Gammaproteobacteria, as described for other herbivorous insects 25 . Gammaproteobacteria has also been described as the dominant class in the only culture-independent study of the microbiome of the genus Monochamus 2 . In the current study, the main bacterial orders detected in the two insect vectors of B. xylophilus were: Enterobacteriales, Pseudomonadales, Vibrionales and Oceanospirilalles. Two OTUs (affiliated to Vibrio and Enterobacteriaceae) were present in all insect samples representing more than 10% of the total number of reads.
Vibrionales and Oceanospirilalles were for the first time reported in association with B. xylophilus insect-vectors. The most represented genera for these orders were Vibrio and Halomonas, respectively. Both Halomonas spp. and Vibrio spp. have been reported to have cellulolytic activity 27,28 . Vibrio species were also detected in the hindgut of Balanogastris kolae, an herbivorous insect pest, with proteolytic and amylolytic activity 29 . This ability for degradation of cellulose, proteins and starch appears to be important for herbivorous insects as the ones in this study.
Both Enterobacteriales and Pseudomonadales orders have already been described in M. galloprovincialis and M. alternatus, in culture-independent and dependent methods respectively 2,20 .
The mentioned bacterial groups may have a role in the insect life cycle, in its different stages, affecting PWD spreading. Egg-and larvae-associated microbiome, which may be vertically transmitted from the adult female, is described as very important for their survival 30 . For instance, Pseudomonas, the most represented genus of Pseudomonadales found in this study, was reported as having a role in insect egg protection against attack by other microbes 31 .
As suggested for B. xylophilus 19 , the insect vector microbiome might also play a role in the defence against harmful compounds produced by the plant. Pseudomonas may help both through supressing plant defences as described by Chung and colleagues (2013) 32 or through chemical detoxification, once this genus was reported to survive the stress of α -pinene (an organic compound found in coniferous trees), being able to use this compound as a carbon source 17 . Members of the same genus were also detected in other studies as part of consortia that protect insects from polluted environments 33 and in other longhorn beetles putatively implicated in the degradation of aromatic hydrocarbons 34 .
Other bacterial groups detected in this study are known to display metabolic capacities that might help the survival of B. xylophilus and its insect vector during pine infection. Vibrio spp. are known to play a role in processing pollutants, including aromatic compounds 35,36 . Stenotrophomonas, represented by OTU 10 present in all sampling locations, was also reported in previous studies to survive the stress of α -pinene being able to degrade this compound 17 .
Besides Monochamus spp., Enterobacteriales and Pseudomonadales were also previously described in B. xylophilus isolated from different pine species (including P. pinaster and P. densiflora) and also in Pinus pinaster trees 2,37-39 . The presence of the same orders in the three organisms that take part in the cycle of the PWD suggests that these orders might have a role in the disease mechanism. Some groups of these bacterial orders are considered quarantine pests (that may represent a risk with potential economic importance) 40 , mostly for their toxin-producing ability. Members of Enterobacteriales were described as producers of cell wall degrading enzymes (cellulases, peptidases and chitinases) and as highly resistant to H 2 O 2 (the most predominant reactive oxygen species, ROS, produced by the plant-basal defences against biological stresses) 41,42 . These abilities may be relevant in the PWD mechanism.
The presence of the same bacterial groups in the three PWD organisms may also be related to the transference of bacteria between these organisms. For example, the origin of B. xylophilus bacterial community is unknown and in previous studies it has been suggested that it might be transmitted from the insect-vector to the nematode 2,23 but more studies must be conducted to confirm this theory.
Both methods used in this study, DGGE and pyrosequencing, revealed differences between the trachea microbiomes of M. alternatus and M. galloprovincialis, suggesting the existence of a species-specific microbiome. Vibrionaceae, Pseudomonadaceae and Halomonadaceae are highly abundant in M. galloprovincialis samples having a minor importance in the communities from M. alternatus. Inversely, the order Xanthomonadaceae is more abundant in M. alternatus samples and the family Acidobacteriaceae is present only in M. alternatus samples, being the third most abundant family in this species. In general, pyrosequencing revealed higher diversity in M. alternatus than in M. galloprovincialis. The existence of a species-specific microbiome was also reported for other insect species 43 . Bacterial communities might have to adapt to differences in morphology and life cycle between species with the consequent variations in the community structure. Also, these two insect species do not co-exist in the same country, having different tree hosts. As reported for B. xylophilus 44 , the bacterial community associated to Monochamus spp. may also differ according with the environmental factors and pine species host. According to the Köppen and Geiger classification 45 , Hikobe (the sampling site in Japan) is Dsa (D -Cold; s -dry summer; a -hot summer) with cold temperate weather and both Portuguese locations are Csa (C -Temperate) with warm temperate weather. Also, different pine species from different locations have differences in their physiology that would affect their microbiota and consequently the insects they host 46 . Thus, these Monochamus species must be adapted to their endemic regions and their tree hosts, with differences in their bacterial community to favour that adaptation.
The hypothesis that the differences between the bacterial communities from the two insect species here studied are due to the unbalanced proportion of insects with and without B. Lactobacillales in M. alternatus bacterial community. Despite the existence of several culture-independent studies describing the bacterial community associated to B. xylophilus 17,47 , members of the orders Acidobacteriales and Lactobacillales were not reported so far in association with the nematode. Thus, the abundance of these orders in the microbiome of M. alternatus does not seem to be related to the presence of B. xylophilus.
Concerning M. galloprovincialis samples, there were no significant differences between the microbiome of insects collected in Comporta and Mortágua, suggesting that sampling location had no influence in M. galloprovincialis microbiome. This stability in the microbiome may suggest relevance in the disease mechanism. Comporta and Mortágua are both temperate regions and the same Pinus species was collected (Pinus pinaster). However, there were site-specific OTUs, like OTU 192 affiliated with Pseudomonas, which was present only in Mortágua samples in high abundance. These specificities may result from the bacterial community adaptation to the local environment; despite the mentioned clime resemblances, sampling locations have differences in altitude (locations differ in 94 m of altitude) and precipitation regimens.

Conclusions
This study suggests that the microbiome of Monochamus tracheae is species-specific and independent of the insect gender and sampling location. More studies will be necessary to determine both the effect of the presence of B. xylophilus in the insect microbiome and the effect of the host tree. Several bacterial groups detected here have been described in other studies in association with the other two PWD organisms (Pinus pinaster trees and B. xylophilus) suggesting that these bacteria might play a role in the disease development. Some taxa might be involved in processes of detoxification of harmful compounds, others might produce toxins, both helping the tree invasion. The insect microbiome characterization is important not only to better understand the role of bacteria in PWD but also for the elaboration of bio-control strategies.

Methods
Sampling. Symptomatic adult trees with symptoms from class IV according with Proença and colleagues (2010) 38 -withered trees with more than 80% of brown needles -were selected in the spring (April and May) of 2013 in locations known to be affected by PWD 48,49 : Pinus pinaster trees were selected in Portugal from Mortágua (40°24′ 17.2″ N 8°13′ 32.7″ W) and Comporta (Herdade da Comporta, 38°22′ 48.67″ N/8°47′ 25.00″ W), central and south regions of Portugal respectively, and Pinus densiflora trees in Japan from Hikobe (Shiwa District, Iwate prefecture, 34°27′ .00″ N/133°46′ .00″ E). Specifically, log sections were sampled, after being inspected for the presence of Monochamus entry-points. Logs were placed in black plastic containers covered at the top with semi-transparent cloth meshes and kept inside a greenhouse. The wood material from each tree was divided according to its origin and kept isolated. All insects used for this study were caught alive after emergence from the tree logs. Those morphologically identified as M. galloprovincialis and M. alternatus were sexed.
As described by Kobayashi and colleagues (1984) 12 , Bursaphelenchus xylophilus concentrates mainly in the vectors' tracheal system. For that reason, the bacterial community associated to insects' tracheae may play a role in the PWD mechanism. Thus, each insect was surface sterilized in a 70% ethanol solution for 1 minute. After rinsing in sterile distilled water, insects were dissected under a binocular scope (Wild Heerbrugg, Switzerland) for trachea extraction.
DNA extraction from insect tracheae, molecular analysis for Monochamus species confirmation and B. xylophilus detection. Tracheae were placed in 40 μl of sterile 1× phosphate-buffered saline (PBS) and total DNA extraction was performed using PowerSoil ® DNA Isolation Kit according to the manufacturer's instructions (MOBIO, CA, USA). Insect morphological identification was confirmed by amplification and sequencing of cytochrome c oxidase I (COI) gene using the primers C1-J-2183a (5′ -CAACAYTTATTTTGATTTTTTGG-3′ ) and TL2-N-3014 (5′ -TCCAATGCACTAATCTGCCATATTA-3′) 50 . The amplification conditions were as follows: an initial denaturation step at 95 °C for 10 min, followed by the three step amplification with 30  Amplicons were purified with the JetQuick ® PCR Product Purification Spin Kit (Genomed, Germany) following the manufacturer's instructions and used as template in the sequencing reactions carried out by the company GATC (Germany) with the primer C1-J-2195 (5′ -TTGATTTTTTGGTCATCCAGAAGT-3′) 50 .
All samples were screened for B. xylophilus by PCR amplification of the species specific MspI satDNA according with the European and Mediterranean Plant Protection Organization 51 instructions using the primers and conditions previously described 52 .

Denaturing Gradient Gel Electrophoresis (DGGE) analysis. A nested PCR technique was applied
to amplify the V3 region of the 16S rRNA gene in order to increase sensitivity. Total DNA extracted as described above was used as substrate. First, the 16S rRNA gene was amplified using the universal primers 27F (5′ -AGAGTTTGATCCTGGCTCAG-3′ ) and 1492R (5′ -GGTTACCTTGTTACGACTT-3′) 53 . The amplification conditions and reaction mixtures consisted of an initial denaturation step (94 °C for 3 min), followed by 30 amplification cycles comprising the following steps: denaturation (94 °C for 1 min), annealing (52 °C for 1 min) and extension (72 °C for 2 min), and a final extension step (72 °C for 10 min). The reaction mixtures were performed as described in Henriques et al. 25 . The second PCR was conducted using 1 μL of the first PCR product as template and the same reagents and volumes except for the primers. The V3 region of bacterial 16S rRNA gene was the target of the second amplification using primers and conditions as described in Henriques et al. 25 . Positive and negative controls were used in both PCRs, replacing the template DNA with DNA from Escherichia coli ATCC 25922 and with sterile dH 2 O respectively. PCR products were applied onto 8% polyacrylamide (37.5:1, acrylamide/bisacrylamide) gel with linear denaturing gradient ranging from 30% to 65% (100% corresponds to 7 M Urea and 40% formamide). Electrophoresis was conducted in a DCode System (Bio-Rad) as described in Henriques et al. 25 .
Gel images were acquired using a Molecular Imager ® Gel Doc TM XR+ System (Bio Rad Laboratories, Hercules, California, USA). Every DGGE gel contained, at its ends, two lanes with a standard of eight bands for internal and external normalization and as an indication of the quality of the analysis 25 . Replication of the PCR and DGGE steps was performed to confirm the consistency of the profiles. DGGE patterns were analysed using Bionumerics Software (Applied Maths, Belgium). Cluster analysis of DGGE profiles was performed using the UPGMA method (group average method) applying Pearson correlation coefficient.
Barcoded 454 Pyrosequencing. Samples were prepared for barcoded 454 pyrosequencing by nested-PCR amplification of the 16S rRNA gene. The first amplification reaction was done as described for DGGE analysis; the second PCR targeted the V3-4 hypervariable region of the 16S rRNA gene, using the forward primer 5′ -ACTCCTACGGGAGGCAG-3′ and the reverse primer 5′ -TACNVRRGTHTCTAATYC -3′ 54  Sequences were processed using both UPARSE 55 and QIIME 56 pipelines on a computer using the Linux© operating system.
In the UPARSE workflow, barcodes were striped and reads were quality filtered to a maximum expected error of 1.0 and trimmed to 400 bp. Sequences were dereplicated (identical reads were merged) as described by Edgar (2015) 55 , generating a file containing a set of unique read sequences each marked with and integer value indicating its abundance. After singletons removal, Operational Taxonomic Units (OTUs) were defined at 97% similarity applying the UPARSE-OTU algorithm (that simultaneously identifies and discards chimeras) to the set of unique read sequences. Chimera check was done a second time running uchime_ref on the OTU sequences generated. Taxonomy assignment was made through QIIME using Uclust as assignment method and Greengenes reference databases. Concerning alpha-diversity, rarefaction plot for the observed species for each sample was made using QIIME scripts; richness index, S, Shannon index of diversity, H 57 and the equitability index, E 58 were calculated using PRIMER software 59  n i is the OTU abundance, S is the number OTUs (used to indicate the number of species) and N is the sum of all reads for a given sample (used as estimates of species abundance) 60 .
Previously to statistical and beta-diversity analysis, sequences were first rarefied to the lowest number of reads in the samples using QIIME script single_rarefaction.py.
All sequences obtained are available in the NCBI platform with the accession numbers SRX1246548, SRX1249956.
Statistical Analysis. Statistical analysis was used to determine if any of the following factors had significant influence (P < 0.05) in the composition of bacterial communities: sex of insects, their species, presence of B. xylophilus and sampling sites. Statistical significance for each factor was evaluated through PERMANOVA based on 999 permutations using PRIMER v6 software 59 . Both DGGE and pyrosequencing dataset were transformed using the log(x + 1) transformation and PERMANOVA was performed on Bray-Curtis distance matrixes constructed from the transformed DGGE table and from the rarefied and transformed OTU abundance table.
Non-metric Multidimensional Scaling (NMDS) was applied to the rarefied and transformed (log(x + 1)) OTU abundance table. This analysis was based on Bray-Curtis distance and computed using R Software (www.r-project. org).
Statistical significance of variance in indexes (S, H and E) was evaluated with Student's t-tests when data had a normal distribution (checked using the Shapiro-Wilk test). For non-normally distributed data, the Mann-Whitney test was applied. Shapiro-Wilk, Student's t and Mann-Whitney tests were computed using R Software (www.r-project.org).