Early changes in bacterial communities in wound tissues of Pinus massoniana after inoculation with Bursaphelenchus xylophilus

Summary – The bacterial communities in the wound tissues of Pinus massoniana were analysed by 16S rDNA amplicon sequencing. The results showed that the bacterial richness and diversity changed remarkably whether the wound was inoculated with pine wood nematode (PWN; Bursaphelenchus xylophilus ) or not after 12 h. However, the predominant bacteria Stenotrophomonas , Burkholderiaceae, Pseudomonas , Serratia and Delftia , introduced by PWN in the wound tissues, changed within 6 h. After 6 h of PWN inoculation, the most abundant genus associating with PWN, Stenotrophomonas , failed to colonise the wound tissues, and the abundance of Delftia decreased, while the other representative bacteria, Burkholderiaceae, Pseudomonas and Serratia , from the PWN were markedly enriched. In addition, our study is the ﬁrst to report the association of Serratia liquefaciens with PWN. Predicted functional analyses using the Tax4Fun tool showed that the alterations in bacterial composition also led to shifts in their functional pathways, especially after 12 h of PWN inoculation. These ﬁndings clariﬁed that the bacteria carried by PWN were responsible for the alterations in bacterial communities in the wound tissues and will shed light on the invasion mechanism of PWN.

Pine wilt disease (PWD) is one of the most destructive diseases of pine trees in the world. The causal agent of the disease was initially confirmed as Bursaphelenchus xylophilus (Steiner & Buhrer, 1934), the pine wood nematode (PWN), by Kiyohara & Tokushige (1971) and for years it was thought to be the only pathogenic agent (Mamiya, 1975(Mamiya, , 1983Nickle et al., 1981;Fukuda et al., 1992). However, some scientists studying the mechanisms of PWD attributed the rapid wilting of pine trees to wilt toxins, and the toxins were probably produced by nematode-associated bacteria (Oku et al., 1979(Oku et al., , 1980Kawazu et al., 1996;Han et al., 2003;Zhao et al., 2003;Tan & Feng, 2004;Guo et al., 2006;Yin et al., 2007). As a result, many experiments have been conducted on these bacteria, including surface bacteria and endobacteria of PWN (Zhao et al., 2000;Wu et al., 2013;Tian et al., 2015), as well as bacteria associated with the beetle vectors (Vicente et al., 2013a;Alves et al., 2016) and hosts of PWN (Proença et al., 2017a), aiming to clarify the origin and exact role of bacteria in PWD. As reviewed by Proença et al. (2017b), the bacteria Pseudomonas, Burkholderia, Serratia, Ewingella and Enterobacter were most commonly reported to be associated with PWD in China, Korea, Portugal and the USA, and other genera, such as Stenotrophomonas, Bacillus and Pantoea, were also carried by the PWN. However, despite the intense research and abundant information about the possible functions of these bacteria, the role of associated bacteria in PWD has mostly been tested in vitro, and when these bacteria come into play in PWD progression remains unclear.
In the wild, the PWN is transmitted from tree to tree by insect vectors, mostly from the genus Monochamus (Coleoptera: Cerambycidae), including Monochamus alternatus in East Asia (Mamiya & Enda, 1972;Lee et al., 1990;Ning et al., 2004), M. carolinensis in North America (Linit et al., 1983) and M. galloprovincialis in Portugal (Naves et al., 2001), through feeding wounds (Linit, 1990) or oviposit wounds (Edwards & Linit, 1992). Afterwards, the PWN moves from the trachea to the tail tip of the vectors and is finally transmitted to the wounds of pine twigs (Aikawai & Tiogashi, 1998), where a new infection cycle begins. Nevertheless, the nematodes hardly migrate or migrate slowly on the wound surface (Tamura, 1984) because most of them are trapped in the sticky resin exuded by epithelial cells on the wound surface (Zhao, 2008) and have to survive harmful secondary metabolites produced by defence mechanisms of the host (Cheng et al., 2013) before successfully invading pine tissues. In this initial stage, some bacteria in the wound may already begin to react against the defence metabolites of the host and help PWN tolerate and overcome the resistance of the host. Therefore, determining what species of bacteria are the pioneers that play a role in this stage is important for illuminating the invasion mechanism of PWN.
To answer the questions mentioned above, the aim of this study was to analyse the changes in bacterial communities in wound tissues caused by inoculation of PWN in the initial stage using 16S rDNA amplicon sequencing, aiming to clarify the bacterial pioneers that would affect the successful invasion of PWN.

EXPERIMENTAL MATERIALS
A virulent PWN strain was initially isolated from naturally infected Pinus massoniana (Fujian, P.R. China). After morphological identification, the nematodes were propagated on Pestalotiopsis sp. (Xie et al., 2017a, b;Sriwati et al., 2008) cultured on PDA medium at 28°C for 7 days. The harvested nematodes were then collected using a Baermann funnel. Large quantities of nematodes were obtained based on repetitive propagation using this method and stored at 4°C until use. Three-year-old P. massoniana seedlings with similar growing conditions (height, 70-80 cm; diam., 1 cm), planted in the Institute of Forest Protection in Fujian Agriculture and Forestry University, China, were used in the inoculation experiment.

INOCULATION AND SAMPLING
The cultured nematodes were rinsed three times with sterile deionised water before use and adjusted to 30 μl nematode suspension (approximately 1000 nematodes) in each 0.5 ml PCR tube; this solution was then directly pipetted onto the artificial wounds of the seedlings. Each seedling had five artificial wounds that were generated by scraping the bark with a razor blade. The longitudinal and axial lengths of the wounds were 3 cm and 0.5 cm, respectively, and the wounds were 5 cm apart from each other. The plants inoculated with the same amount of sterile water were used as controls. To distinguish the origin of bacteria in the wound tissues after inoculation, the bacterial communities carried by the PWN and the communities that originally existed in the cortex tissues of healthy P. massoniana seedlings were detected.
As most of the nematodes were trapped on the wound surface for approximately 6-12 h according to different host and inoculum densities (Jin, 2007;Zhang et al., 2007;Li & Ye, 2008;Su et al., 2008), we chose 6 h and 12 h as the experimental time points. The tissues of the wounds, mostly the cortex, were collected 6 and 12 h after inoculation and then stored at −80°C until use, as were the corresponding two control groups. The cortex tissues of healthy P. massoniana seedlings were collected immediately after the bark had been scraped and the rest of the cultured nematodes (three tubes with 0.5 ml nematode precipitate in each) were rinsed three times with sterile deionised water before storing at −80°C. Because of the limited amount of cortex tissues from one single wound, a mixture of 15 cortex tissues from 15 wounds (three seedlings) was treated as a treatment sample, and each treatment had three replicates (nine seedlings). Thus, a total of 45 pine seedlings were used in this research, with nine seedlings in each group (two experimental groups, two control groups and one group of healthy pine tissues). The tools used above were all decontaminated before use, and the experiment was conducted at 28°C in a glasshouse.

DNA EXTRACTION
Total genomic DNA of the samples was extracted using the CTAB method (Lutz et al., 2011), and the concentration and purity were detected by 1% agarose gel electrophoresis. Afterwards, the DNA products were diluted to 1 ng μl −1 with sterile water before PCR amplification.

PCR AMPLIFICATION
To study the diversity and composition of bacteria in the wound tissues, the distinct V3-V4 regions of 16S rDNA were PCR amplified with specific barcoded primers: 341F (5 -CCTAYGGGRBGCASCAG-3 ) and 806R (5 -GGACTACNNGGGTATCTAAT-3 ). PCR was performed in a 30 μl volume, containing 15 μl Phusion ® High-Fidelity PCR Master Mix (New England Biolabs), 0.1 μM each primer and 10 ng template DNA. Thermal cycling consisted of initial denaturation at 98°C for 1 min, followed by 30 cycles of denaturation at 98°C for 10 s, annealing at 55°C for 30 s, 72°C for 30 s and a final extension at 72°C for 5 min.

GENE LIBRARY PREPARATION AND SEQUENCING
PCR products were detected using 2% agarose gel electrophoresis, and then purified with a GeneJET Gel Extraction Kit (Thermo Scientific). Sequencing libraries were generated using the Ion Plus Fragment Library Kit 48 rxns (Thermo Scientific) following the manufacturer's recommendations. After assessing the library quality on the Qubit@ 2.0 Fluorometer (Thermo Scientific), the PCR products were sequenced on an Ion S5™ XL platform, and 600 bp single-end reads were generated.

BIOCHEMICAL ANALYSIS
Raw reads were obtained after single-end reads were assigned to samples based on their unique barcode and truncated by removing the barcode and primer sequences. Quality filtering of the raw reads was performed under specific filtering conditions to obtain high-quality clean reads according to the Cutadapt (Version 1.9.1, http:// cutadapt.readthedocs.io/en/stable/) (Martin, 2011) quality control process. Subsequently, clean reads were finally obtained after comparison with the Silva database (https://www.arb-silva.de/) (Quast et al., 2012) using the UCHIME algorithm (http://www.drive5.com/usearch/ manual/uchime_algo.html) (Edgar et al., 2011) to detect and remove the chimeric sequences.
Sequence analysis was performed by UPARSE software (Version 7.0.1001, http://drive5.com/uparse/) (Edgar, 2013). Sequences with 97% similarity were clustered to the same operational taxonomic units (OTUs). A representative sequence for each OTU was screened for further annotation, and taxonomic information was annotated by means of the Silva database based on the MOTHUR algorithm.
To assess microbial diversity with or without infection by PWN, alpha diversity was calculated, including observed species, Chao 1, Shannon, Simpson and good coverage, all of which were calculated with QIIME (Version 1.7.0) (Caporaso et al., 2010) and displayed with R software (Version 2.15.3).
To evaluate the differences in species complexity among samples, beta diversity and unweighted UniFrac diversity were analysed by QIIME software (Version 1.7.0) (Caporaso et al., 2010). The significant differences in bacterial community structure among samples were analysed using MetaStat (White et al., 2009) andLEfSe (Segata et al., 2011). The functional diversity of bacteria in each group was predicted by Tax4Fun (Aßhauer et al., 2015).

BACTERIAL DIVERSITY ANALYSIS
To assess the possible alterations in bacterial diversity caused by PWN inoculation, alpha diversity was analysed based on observed species, Chao 1 (richness), Shannon index and Simpson index (diversity) ( Table 1). All indices indicated that the bacterial richness and diversity in Group Bx6h were comparable with those of its control group (C6h), while Group Bx12h had significantly lower richness and diversity indices compared to the control group (C12h) (P < 0.05, Wilcox). Even in the two control groups, the bacterial richness and diversity changed, with the longer time group having greatly increased indices. Based on the sequencing data, we found that the main genera that increased in abundance (Group C12h compared to Group C6h) was Rhizobacter and other unidentified genera. Moreover, compared to all the other groups, Group Bx had the lowest bacterial richness and diversity (P < 0.05, Wilcox).

BACTERIAL COMMUNITY COMPOSITION
To gain insights into the differences in bacterial communities, we further analysed the individual bacterial taxa in each group and presented the results according to the clusters of Unweighted Pair-group Method with Arithmetic Means (UPGMA) Clustering. As shown in Figure 2, the PWN-inoculated groups (Bx6h, Bx12h) shared a similar relative abundance of bacterial phyla and clustered together, while the bacterial phyla of Group Bx that we used for inoculation were distinctive from the other five groups.

GROUPS TO CONTROL GROUPS
In Group Bx6h, the most represented genus was Methylobacterium (15.5 ± 4.29%), followed by Pseudomonas (8.0 ± 3.94%). With the time of inoculation prolonged (Group Bx12h), the structure and relative abundance of bacteria at the genus level showed no significant difference, except for the remarkable decrease in Nocardioides, which was detected in Group Pm but not in Group Bx.
To identify the change in bacterial communities after inoculation with PWN on the wound surface of P. massoniana, we compared the inoculation groups with their corresponding control groups and generated a heatmap together with relative abundance histogram at the genus level (Fig. 4). The comparative results between Groups Bx6h and C6h showed that the most abundant genus from PWN, Stenotrophomonas, completely failed to colonise the wound tissues, and the abundance of the other rep-resentative genus, Delftia, also dramatically decreased. Nevertheless, the other top represented genera of PWN (unidentified Burkholderiaceae, Pseudomonas and Serratia) were significantly enriched (P < 0.05, MetaStat), accounting for 2.3 ± 0.99%, 8.0 ± 3.95% and 0.9 ± 0.45% of the total reads in this inoculation group, respectively. Comparative analysis between Group Bx12h and the corresponding control group led to similar results except for the remarkably lower relative abundance of Nocardioides.
LEfSe (LDA Effect Size), a software for discovering significant differences between groups, was also used, and the results indicated that both PWN-inoculation groups exhibited elevated proportions of Pseudomonas, unidentified Burkholderiaceae and Serratia (particularly Serratia liquefaciens) (Fig. 5). In addition, although with low relative abundance, other genera, such as Mesorhizobium (belonging to the family Burkholderiaceae), Rhizobacter (controversial taxonomy; our sequencing result classified it as Burkholderiaceae) and Friedmanniella, were also enriched in Group Bx6h, while Yersinia sp. Ha77 and Ewingella americana, both of which belong to Enterobacteriaceae, were enriched in Group Bx12h. Fig. 4. Heatmap and relative abundance of genera (more than 0.5% of the total reads) with significant differences (P < 0.05; MetaStat). The red and black font indicate the significant increase and decrease in abundance of the genera in Group Bx6h and Bx12h, respectively, compared to their control groups, and the up and down arrows in the heatmap demonstrate the specific trend in each pine wood nematode (PWN; Bursaphelenchus xylophilus)-inoculation group. The histogram shows the average relative abundance of these genera in each group (n = 3). Error bars represent SE of the mean. Bx6h, the 6 h PWN-inoculation group; Bx12h, the 12 h PWN-inoculation group; C6h, control group of Bx6h; C12h, control group of Bx12h; Pm, healthy cortex of Pinus massoniana; Bx, B. xylophilus.

FUNCTIONAL PREDICTION
We used Tax4Fun based on the 16S Silva database to predict the functional diversity of the bacteria in each group. A total of 6471 functional orthologues were predicted, assigning to 43 level 2 KEGG pathways, with membrane transport, carbohydrate metabolism, amino acid metabolism and translation having greater abundance of related genes. A heatmap was made from the top 25 KEGG orthologue (KO) groups according to functional annotation and abundance, and all groups were clustered based on functional relative abundance. As expected, the functions of the bacteria associated with Group Bx were significantly different from those associated with the other groups (data not shown), with higher KO abundance in metabolism-related pathways such as 'enzyme families' and 'glycan biosynthesis and metabolism', environmental information processing pathways such as 'membrane transport' and 'signal transduction', and cellular processes pathways such as 'cell motility' and 'cellular community prokaryotes'. Among the other groups from the wound tissues of P. massoniana, the bacterial functions of Group Bx6h were more similar to those of its control group and the other two non-PWN-inoculation groups, while the bacterial functions of Group Bx12h clustered farther from all the other groups (Fig. 6). As seen from the heatmap, Group Bx12h had higher KO abundance in 'carbohydrate metabolism', 'xenobiotic biodegrada-tion and metabolism', 'amino acid metabolism', 'lipid metabolism', 'metabolism of terpenoids and polyketides', 'nucleotide metabolism', 'enzyme families', 'transport and catabolism' and 'replication and repair' pathways. In level 2 KEGG pathways, both Group Bx6h and Group Bx12h had no significant alteration in KO abundance compared to their control groups. However, Group Bx12h had significantly higher KO abundance in the 'metabolism of cofactors and vitamins' and 'nucleotide metabolism' pathways than Group Bx6h.

Discussion
Generally, the bacterial populations and diversity inside trees increase with disease progression , but little is known about the bacteria in wound tissues in the initial stage before the pathological response of pines to PWN invasion can be seen. Our study illustrates the changes in bacterial communities in the wound tissues of P. massoniana after inoculation with PWN. The results showed that the bacterial richness and diversity changed little after 6 h but decreased significantly after 12 h compared to their control groups. However, obvious alterations in some of the predominant bacteria carried by PWN occurred within 6 h; in particular the unidentified Burkholderiaceae, Pseudomonas and Serratia ge- nera were enriched but the abundance of Delftia and Stenotrophomonas were either reduced or excluded.
Unlike the increasing population dynamics of bacteria in disease stages, the bacterial colonies in the adjacent xylem of inoculation sites hardly changed after 3 h and were generally lower than those in the samples with disease symptoms (Roriz et al., 2011), probably because there is a barrier outside the xylem, which is the defence response of the host, leading to little change in bacterial communities in the xylem. In our study, the inoculation of PWN apparently directly changed the bacterial communities in the wound tissues of PWN, particularly after 12 h, which may be enough time for adaptive bacteria to survive and proliferate and for host pine to eliminate unsuitable bacteria. Interestingly, Group C12h had the highest bacterial richness and diversity, indicating that even without PWN inoculation, the bacterial communities in the wound tissues might change. Although the majority of newly added bacteria remained unidentified in this group, the increasing richness was partly due to Rhizobacter, which was also observed in healthy Group Pm but with lower abundance. In addition, Rhizobacter was more abundant in both PWN-inoculation groups than Group Pm; thus, we suspect it might be beneficial for pine trees to confront injuries or resist the invasion of PWN. In support of this suggestion, Dhanasekar & Dhandapani (2012) proposed that Rhizobacter is a nitrogen-fixing bacteria and is mainly involved in the biological control of pathogens, nutrient cycling and seedling establishment.
In Group Bx, the dominant bacterial genera were Stenotrophomonas, unidentified Burkholderiaceae, Pseudomonas, Serratia, Delftia, Novosphingobium and Chryseobacterium. Among these bacteria, Stenotrophomonas was the most dominant genus, in accordance with previous studies (Tian et al., 2010;Wu et al., 2013), and Stenotrophomonas was also detected in PWN obtained from P. massoniana. Pseudomonas, Serratia and Burkholderiaceae have been widely reported in association with PWN and are able to induce PWD symptoms (Oku et al., 1980;Han et al., 2003;Guo et al., 2006Guo et al., , 2007Proença et al., 2010;Vicente et al., 2011Vicente et al., , 2012. Meanwhile, Serratia, namely, S. marcescens, was associated with M. alternatus, and secondary metabolites of S. marcescens are able to degrade lignin (Fu, 2017). In our study, we also detected S. marcescens in PWN, but with a low proportion (0.1 ± 0.02%). Cheng et al. (2013) and Wang et al. (2019) found that Delftia and Novosphingobium were dominant in PWN isolated in Zhejiang Province in China, and Liu et al. (2017) identified Delftia, Pseudomonas, Stenotrophomonas and Rhizobium as the dominant species on the surface of PWN in the USA. Chryseobacterium was also reported to be one of the main species carried by the PWN from Japan (Ju et al., 2008) and the USA (Proença et al., 2014). It is plausible that bacteria carried by PWN varied in different studies because of multiple factors, such as geographic area, Pinus species, feeding fungus, isolation methodologies and surface sterilisation. To avoid eliminating the potential PWD-related bacteria, the PWN we used here were obtained without sterilisation because a considerable number of bacterial species are present on the cuticle surface of PWN (Roriz et al., 2011;Vicente et al., 2011), and superficially associated bacteria can also potentially be involved in PWD (Li, 2008).
In Group Pm, the bacterial community was dominated by Proteobacteria and Actinobacteria at the phylum level. The top three genera were Methylobacterium, Sphingomonas and unidentified Beijerinckiaceae, all of which belong to Proteobacteria. A recent study also reported that the most dominant endophytic bacteria in healthy, as well as in diseased, P. massoniana is Proteobacteria (Li et al., 2018). However, the corresponding genera were different, and the second most dominant phylum was Bacteroides, which ranked third in our study. These differences could be explained by differences in sampling location, sampling area and tree age. In addition, several studies have confirmed the existence of endophytic bacteria in pine trees, with Methylobacterium being the major genus (Pirttilä et al., 2000(Pirttilä et al., , 2008, possibly simply because it is a beneficial endophytic bacteria for the growth of pine seedlings (Pohjanen et al., 2014).
After inoculation with PWN, unidentified Burkholderiaceae, Pseudomonas and Serratia (especially S. liquefaciens) were introduced to the wound surface together with PWN and significantly enriched within at least 6 h, while interestingly the most dominant genus, Stenotrophomonas, could not colonise wound tissues and the abundance of one of the representative genera, Delftia, was reduced, indicating that unidentified Burkholderiaceae, Pseudomonas and Serratia adapted to the environment containing secondary metabolites from the host but Stenotrophomonas and Delftia did not. For the above result, we may doubt that the increases in Burkholderiaceae, Pseudomonas and Serratia resulted from the in-troduction of PWN, but the abundance of these bacteria also increased in Group Bx12h compared to Group Bx6h, confirming that they were indeed enriched in the wound tissues and may play a role in the initial infection activities. In summary, these bacteria could be considered as biomarkers to understand the development of PWD in this specific stage.
As terpenoids (especially α-pinene) are defence compounds of pine trees and are detrimental to the reproduction of the PWN (Kong et al., 2007), α-pinene degradation is important for the survival and invasion of PWN. It was reported that Pseudomonas is one of the main strains capable of degrading α-pinene; it not only survives the stress of α-pinene but also utilises α-pinene and other secondary metabolites of pine, including benzoate, as a carbon source for growth, whilst Stenotrophomonas cannot survive under the stress of benzoate and grows slowly in LB medium containing α-pinene (Cheng et al., 2013). A recent study observed similar results; the abundance of Pseudomonas sp. carried by PWN increased, while the abundance of Stenotrophomonas and Delftia spp. decreased when PWN was fumigated with α-pinene (Wang et al., 2019), which may explain the enrichment of Pseudomonas and the reduction or even exclusion of Delftia and Stenotrophomonas in our study. However, we could not ignore the fact that one species of Stenotrophomonas, namely, S. maltophilia, was related to the strong virulence of PWN (Wu et al., 2013), and when S. maltophilia NSPmBx03 was treated with PWN, many cell wall degradation-related and detoxificationrelated genes, including pectate lyase, glutathione Stransferase, ATP-binding cassette transporter and cytochrome P450, might be upregulated, and the virulence of PWN might be enhanced (He et al., 2016). In addition, some Stenotrophomonas spp. were reported as being able to degrade aromatic compounds and other xenobiotics (Mangwani et al., 2014;Tiwari et al., 2016), which seems beneficial for PWN to invade its host. However, other strains of S. maltophilia might inhibit PWN hatching (Tian et al., 2015) or possess nematicidal activity against PWN (Huang et al., 2009), leading to the opposite result. A similar controversy was also presented by Serratia: many species exhibited high toxicity against PWN in vitro (Paiva et al., 2013) but some were able to produce copious cellulose, form biofilms, resist reactive oxygen species, tolerate growth in the presence of xenobiotic/organic compounds and not only colonise themselves (Nascimento et al., 2018;Vicente et al., 2012Vicente et al., , 2013aVicente et al., , 2016b but also assist PWN survival under prolonged oxidative stress con-ditions (Vicente et al., 2013b(Vicente et al., , 2016b. Considering these contradictory results on different strains within the same genera, the role of these bacteria in PWN and PWD needs more investigation. Interestingly, within the Serratia genus, S. liquefaciens was reported for the first time in our study to be related to PWN. It was relatively abundant (3.3%) in Group Bx and enriched in the wound tissues after 6 h of PWN inoculation. Although its detailed role in PWN remains to be further studied, this species has been described to produce volatile organic compounds and thus could promote plant growth and attract Caenorhabditis elegans (Grewal & Wright, 1992;Wang, 2018).
To predict the functions of bacteria, we used the Tax4Fun approach. Because of the limitation of this method, we could not identify the exact role of each bacterium, but the results of this section could suggest that the alterations in bacterial composition potentially lead to shifts in their functional pathways. This result was more obvious in the longer inoculation time group (Group Bx12h) in our study and was reflected mainly in the increase in KO abundance in multiple metabolism-related pathways such as 'carbohydrate metabolism', 'xenobiotics biodegradation and metabolism', 'amino acid metabolism', 'lipid metabolism', and 'metabolism of terpenoids and polyketides', which potentially indicates that related bacteria may play a role in the invasion of PWN at 12 h. Further studies are needed to determine the specific role of these bacteria in pine trees in vivo.