Bacteriome from Pinus arizonica and P. durangensis: Diversity, Comparison of Assemblages, and Overlapping Degree with the Gut Bacterial Community of a Bark Beetle That Kills Pines.

Symbioses between plants and microorganims have been fundamental in the evolution of both groups. The endophytic bacteria associated with conifers have been poorly studied in terms of diversity, ecology, and function. Coniferous trees of the genera Larix, Pseudotsugae, Picea and mainly Pinus, are hosts of many insects, including bark beetles and especially the Dendroctonus species. These insects colonize and kill these trees during their life cycle. Several bacteria detected in the gut and cuticle of these insects have been identified as endophytes in conifers. In this study, we characterized and compared the endophytic bacterial diversity in roots, phloem and bark of non-attacked saplings of Pinus arizonica and P. durangensis using 16S rRNA gene pyrosequencing. In addition, we evaluated the degree of taxonomic relatedness, and the association of metabolic function profiles of communities of endophytic bacteria and previously reported gut bacterial communities of D. rhizophagus; a specialized bark beetle that colonizes and kills saplings of these pine species. Our results showed that both pine species share a similar endophytic community. A total of seven bacterial phyla, 14 classes, 26 orders, 43 families, and 51 genera were identified. Enterobacteriaceae was the most abundant family across all samples, followed by Acetobacteraceae and Acidobacteriaceae, which agree with previous studies performed in other pines and conifers. Endophytic communities and that of the insect gut were significantly different, however, the taxonomic relatedness of certain bacterial genera of pines and insect assemblages suggested that some bacteria from pine tissues might be the same as those in the insect gut. Lastly, the metabolic profile using PICRUSt showed there to be a positive association between communities of both pines and insect gut. This study represents the baseline into the knowledge of the endophytic bacterial communities of two of the major hosts affected by D. rhizophagus.

Introduction

Symbiosis has been recognized as a “key driver” force in the species evolutionary process. Plant-microorganism interactions have resulted in beneficial, neutral or detrimental effects for both groups. In particular, bacterial communities carry out different biological activities in organs and tissues of plants: rhizosphere, phyllosphere, and endosphere (Turner et al., 2013). Endophytic microorganisms of the endosphere, i.e., those that reside within the inner tissues of plants without causing any symptoms of disease (Hallmann et al., 1997), have enhanced plant tolerance to different biotic and abiotic factors and mediated the interaction with parasites (Reinhold-Hurek and Hurek, 2011; Hardoim et al., 2015; Hashem et al., 2016).

Research on microbial endophytes in plants has been mainly focused on crops of economic importance (e.g., corn, sorghum, soybean, wheat) (Zinniel et al., 2002), but little attention has been paid to wild plants, especially conifers (Carrell and Frank, 2014). A number of studies in conifers have partially characterized the endophytic diversity in different tissues (roots, stems, buds, needles, seeds, and pollen), but almost none have explored the ecological and functional role of these bacteria (O’Neill et al., 1992; Shishido et al., 1995; Pirttilä et al., 2000; Strzelczyk and Li, 2000; Cankar et al., 2005; Izumi et al., 2008; Bal et al., 2012). In addition, most of these studies have been carried out using culture-dependent methods and traditional molecular techniques (e.g., molecular cloning, fingerprinting). Given that these approaches are limited in their statistical coverage, they can lead to biases in diversity characterization by the low number of sequences recovered. In contrast, the application of next generation sequencing (NGS) technologies overcomes these limitations and allows a better taxonomic and functional characterization of endophytic communities in conifers (Carrell and Frank, 2014, 2015; Carrell et al., 2016; Kaul et al., 2016; Rúa et al., 2016).

The Larix, Pseudotsugae, Picea, and Pinus genera are hosts of a group of weevils, Dendroctonus bark beetles, because they carry out their life cycle inside trees. These insects play an important role in coniferous forests as agents of natural renovation and regeneration, but they are also one of the most destructive insects of these communities because they kill from 100 to 1000s of healthy trees during periodic outbreaks (Wood, 1982; Reeve et al., 2012). The life cycle of these bark beetles is apparently simple. Females locate and arrive on trees following host kairomones (Gray et al., 2015). Later, they bore the bark toward the phloem, and once inside produce pheromones to attract males, sometimes leading to mass attacks on trees. Pairs copulate and excavate galleries where females oviposite. Then, larvae hatch, feed on the stem phloem, develop into new adults and lastly, emerge to colonize other pine trees (Raffa et al., 2008).

Several studies have documented the diversity and functional role of bacterial communities associated with certain Dendroctonus species and their galleries using different methodologies (Cardoza et al., 2006; Scott et al., 2008; Adams et al., 2009, 2013; Morales-Jiménez et al., 2009, 2012, 2013; Boone et al., 2013; Durand et al., 2015; Mason et al., 2015; Dohet et al., 2016; Xu L.T. et al., 2016; Xu L. et al., 2016; Briones-Roblero et al., 2017a,b). In addition, Hernández-García et al. (2017) determined the presence of a strict core bacteriome in the gut of these insects, constituted by Enterobacter, Pantoea, Pseudomonas, Rahnella, Raoultella, and Serratia, and a relaxed core composed of Acinetobacter, Propionibacterium, Providencia, Stenotrophomonas, Erwinia, Kluyvera, Paenibacillus, and Ralstonia. Several of these bacterial genera present in the gut and cuticle of Dendroctonus bark beetles have been reported as endophytes in pines and other coniferous genera (Pirttilä et al., 2000; Strzelczyk and Li, 2000; Cankar et al., 2005; Izumi et al., 2008; Bal et al., 2012; Carrell and Frank, 2015; Rúa et al., 2016), suggesting that these bacteria might be environmentally acquired, mainly during phloem feeding (Hernández-García et al., 2017).

In this study, we characterized and compared the endophytic bacterial diversity in roots, phloem and bark of non-attacked saplings of Pinus arizonica Engelm (Arizona pine) and P. durangensis Martinez (Durango pine), using 16S rRNA gene pyrosequencing. In addition, we explored the degree of taxonomic relatedness and the association of functional metabolic profile of endophytic bacterial communities and previously reported gut bacterial communities of the bark beetle, D. rhizophagus (Briones-Roblero et al., 2017a). These pine species were selected based on the fact that they are two preferential hosts of this insect which colonizes and kills saplings of 11 pine species, mainly Arizona and Durango pines, in Mexico (Salinas-Moreno et al., 2004; Mendoza et al., 2011). This is also because the gut bacterial community of this insect has been widely characterized and a core bacteriome has been determined from all life stages (Briones-Roblero et al., 2017a).

Materials and Methods

Sample Collection and Surface Sterilization

Healthy saplings of Arizona and Durango pines (<3 m tall, 10 cm diameter), non-attacked by Dendroctonus rhizophagus, were collected in San Juanito, Bocoyna Municipality, Chihuahua state (27° 45′ 11′′ N 107° 38′ 06′′ W, 2288 masl), and El Salto locality, Pueblo Nuevo Municipality, Durango state, Mexico (23° 41′ 31′′ N 105° 43′ 19′′ W, 2702 masl), respectively. Six trees of each pine species were removed from the soil using a peak and shovel, to integrate two biological replicates of three trees each one. Roots and stems were covered with plastic film, and stored at 4°C for transportation to the laboratory where they were immediately processed.

Plant parts were surface rinsed with sterile water for eliminating rhizospheric soil and contaminants. Several small fragments of approximately 2 cm × 3 cm of roots, phloem, and bark were obtained with sterile knife and fine forceps, placed in sterile resealable bags and stored at -20°C until processing. Surface sterilization for eliminating epiphytic microorganisms was performed according to the methods described by Mendes et al. (2007) with modifications to timing and concentrations. Briefly, fragments of each tissue were immersed separately in 70% ethanol for 3 min, 4% NaClO for 5 min, and 70% ethanol for 30 s; and finally, samples were rinsed in deionized water six times. The surface sterilization of tissues was confirmed by no growth of bacteria in Petri dishes with tryptic soy agar (TSA, BD, Difco, United States) inoculated with the last rinsing water and by negative PCR amplification of that same water.

Metagenomic DNA Extraction, 16S rRNA Gene Amplification, and Pyrosequencing

We extracted DNA from 36 samples including three samples of root, phloem, and bark of each biological replicate from each pine species (Arizona and Durango pines). Tissues were independently ground to fine powder with sterilized mortars, pistils, and liquid nitrogen in a sterile environment. DNA of each sample was extracted using the protocol described by Zhou et al. (1996) with minor modifications. Briefly, 1 mL of lysis buffer (100 mM Tris-HCl pH 8, 100 mM EDTA pH 8, 100 mM sodium phosphate pH 8, 1.5 M NaCl, 1% CTAB, 1.2% Triton X-100) was added to ∼50 mg of each sample contained in 2 mL screw-cap tubes. Suspensions were vortexed for 10 min and 30 μL of lysozyme was added, after which samples were homogenized for 20 s and incubated at 37°C for 1 h. Proteins were degraded by the addition of 30 μL of proteinase K (10 mg/mL), vortexing for 10 s, and incubating at 60°C for 4 h with gentle end-over-end inversions every 30 min. The aqueous phase was collected after centrifugation at 14,000 × g for 2 min at room temperature (RT) into a new 2 mL sterile tube. Proteins were removed by the addition of an equal volume of phenol-chloroform-isoamyl alcohol (25:24:1, v/v/v) and the aqueous phase was recovered by centrifugation at 10,000 × g for 10 min (this step was repeated twice). The DNA was precipitated with the addition of an equal volume of cold isopropanol, overnight incubation at -20°C, and centrifugation at 10,000 × g for 20 min. The DNA pellet was washed twice with 250 μL cold ethanol at 10,000 × g for 5 min, and the ethanol was then discarded and the tube dried at RT, finally being resuspended in sterile deionized water. DNA was quantified using a NanoDrop 2000c Spectrophotometer (Thermo Scientific, Wilmington, DE, United States) and was observed in 1.0% agarose gel.

DNA from 36 samples was amplified independently with 8F and 556R primers (Navarro-Noya et al., 2013), tagged with 10 bp multiplex identifiers barcode (MID) and a Roche 454 adaptor (Roche, Mannheim, DE, United States) according to the Lib-L protocol. The V1-V3 region of the 16S rRNA gene was amplified as reported by Briones-Roblero et al. (2017a). Amplification products were purified using a QIAquick Gel Extraction kit (Qiagen, Valencia, CA, United States) according to the manufacturer’s protocol and quantified using a Nanodrop 2000c spectrophotometer. For pyrosequencing using a Roche GS-FLX Titanium 454 pyrosequencer (Roche, Mannheim, DE, United States) at Macrogen, Inc., (Seoul, South Korea), we pooled three amplification products from each tissue in equimolar concentrations (100 ng/μL), for a total of 12 libraries (Table ).

Analysis of Pyrosequenced Data

The raw sequencing data were analyzed using the Quantitative Insights Into Microbial Ecology (QIIME) pipeline v.1.8 (Caporaso et al., 2010). Low-quality reads were eliminated according to the follow filtration criteria: Phred quality score < 25, sequences < 150 or > 500 bp long, sequences that had homopolymers > 6, sequences that had any ambiguous characters, and any errors in barcodes.

High-quality sequences were sorted into operational taxonomic units (OTUs) at 97% of similarity threshold using the UCLUST algorithm (Edgar, 2010) based on the open-reference method. Chimeric sequences were detected and eliminated using Chimera Slayer (Haas et al., 2011). An OTU table was constructed with representative sequences for each OTU (longest) and their relative abundance across the samples. Sequences matching chloroplasts were manually removed. All representative sequences were aligned using PyNast (Caporaso et al., 2010) and the taxonomic assignment to the different levels, from phylum to genus, was performed at an 80% confidence threshold using the Ribosomal Database Project (RDP) Naïve Bayesian Classifier[1] (Wang et al., 2007). In order to corroborate the taxonomic assignment of OTUs, they were manually compared against reference sequences deposited in three databases: RDP[2], GenBank[3] (Altschul et al., 1990), and Greengenes[4] (DeSantis et al., 2006).

Samples were homogenized with respect to the sample with the lowest number of reads. The α- and β diversities were estimated in QIIME for each sample (Kuczynski et al., 2011). For estimation of α diversity, the replicates of each tissue of each pine species were averaged and used the following metrics: Chao1 (richness estimator), Shannon, Simpson, and Phylogenetic Diversity (PD) (diversity indices) (Magurran, 1998; Faith and Baker, 2007). Rarefaction curves and the Good’s coverage index were computed to determine sampling completeness (Chao et al., 1988). Owing to unequal variances, differences in the means of richness and diversity values were tested using Welch’s test (ANOVA).

A heatmap was constructed with the relative abundance information of bacterial genera identified and merged with a dendrogram built through the unweighted pair group method with arithmetic mean (UPGMA) method using the Bray–Curtis similarity index. The clusters reliability was supported by bootstrap test after 1000 pseudoreplicates and by the cophenetic correlation coefficient estimated in PAST v.3.11 (Hammer et al., 2001); the tree was visualized in iTol[5] (Letunic and Bork, 2007).

Based on presence/absence data of bacterial genera, Venn diagrams were generated with the web application, Venny (Oliveros, 2007), for the following comparisons: (a) between the total community of pine species; (b) among tissue samples of each pine species; and (c) among the total community of each pine species and that previously reported for D. rhizophagus (Briones-Roblero et al., 2017a).

The β diversity of endophytic bacterial communities of each pine species was estimated using the Bray–Curtis dissimilarity index as well as unweighted and weighted Fast UniFrac distances metrics. Both UniFrac estimators are based on phylogenetic richness, but weighted UniFrac analysis also considers the relative abundance of OTUs (Lozupone et al., 2011). This information is extracted from a phylogenetic tree constructed with the maximum likelihood method in FastTree (Price et al., 2009), employing the generalized time-reversible (GTR) nucleotide model. A Monte Carlo test was applied to determine statistically significant differences in β-diversity values among bacterial communities after 1000 randomizations of the original UniFrac matrices. To compare endophytic bacterial communities among different tissues of Arizona and Durango pines, along with the gut community of D. rhizophagus (Briones-Roblero et al., 2017a), we carried out principal coordinate analyses (PCoA) using the Bray–Curtis dissimilarity index and both UniFrac (unweighted/weighted) indexes; the corresponding plots of these analyses were visualized in NTSYSpc v.2.02 (Rohlf, 1997). Analyses of similarity (ANOSIM) were conducted to establish significant differences in the community composition among the analyzed tissues, as well as between endophytic and gut bacterial communities.

To determine the degree of taxonomic relatedness (i.e., the relative closeness of species in the taxonomic hierarchy) of the most abundant members found in the gut of D. rhizophagus (Acinetobacter, Rahnella, Serratia, Stenotrophomonas, Propionibacterium, and Pseudomonas) and their equivalents in endophytic bacterial communities, we performed phylogenetic inference analyses of these bacteria. Sequences of bacterial communities from both pines and D. rhizophagus gut, along with certain reference sequences of the GenBank, including some of isolates of this insect, were aligned in Clustal X v.2.1 (Larkin et al., 2007) and trimmed at the 5′ and 3′ ends using SeaView v.4.4.2 (Gouy et al., 2010). The alignments were submitted to the PhyML Web server[6] to select the most appropriate models of nucleotide evolution (Lefort et al., 2017) and perform maximum-likelihood phylogenetic analyses (Guindon et al., 2010). Bootstrap tests were conducted to determine clustering strength and trees generated were rooted with Anabaena variabilis as the outgroup. Lastly, trees were visualized in iTol5 (Letunic and Bork, 2007).

Predictive Functional Profiling

The Phylogenetic Investigation of Communities by Reconstruction of Unobserved States (PICRUSt) software (Langille et al., 2013) was utilized through the web application Galaxy[7] (Afgan et al., 2016) to investigate the predictive functional profile of bacterial communities of different tissues of each pine species. The sequences from each library were demultiplexed, and OTUs were clustered against the Greengenes database with the closed-reference method employing QIIME to generate an OTU table in biom-format. The accuracy of metagenome predictions was determined with the weighted nearest-sequenced taxon index (weighted NSTI). The NSTI summarizes the extent to which microorganisms in a sample are related to sequences genomes and represents the average branch length that separates each OTU in a sample from a reference bacterial genome, weighting their relative abundance in each sample. Low values of this index indicate a closer mean relationship. The table containing the predicted gene family-counts per sample, based on orthologous groups and identifiers following Kyoto Encyclopedia of Genes and Genomes (KEGG)[8] at levels 2 and 3, was cleaned according to these criteria: removal of categories unrelated with bacterial physiology/metabolism (like human diseases) and removal of gene family categories with a count equal to 0. The database at level 2 was used to generate a heatmap in CIMminer[9] and it was merged with a dendrogram built through UPGMA method using the Gower similarity index. Tree reliability was supported by a bootstrap test after 1000 pseudoreplicates and by the cophenetic correlation coefficient estimated in PAST v.3.11 (Hammer et al., 2001). The database at level 3 was employed to build a bar chart with the relative abundance information of the metabolic pathways. Spearman’s correlation was performed to associate the abundances of inferred metabolic functions of pine endophytes and D. rhizophagus gut communities obtained in PICRUSt, at the level 2 of metabolism category, with PAST v.3.11 (Hammer et al., 2001).

Data Accessibility

The Standard Flowgram Format (SFF) file was deposited in the NCBI SRA database under bioproject PRJNA395769 with biosample SAMN07414558 and accession number SRR5872379.

Results

Pyrosequencing Data

A total of 529 bacterial OTUs were identified from 23,330 high-quality and non-plastid sequences in 11 of 12 samples. A biological replicate corresponding to the bark of Arizona pine was discarded given that it presented an insufficient number of reads. The remaining 11 samples featured an average length of 335 bp/read, and an average of 2,120 reads/sample (Table ). The bark of both pine species exhibited the highest OTUs number, and meanwhile, phloem had the lowest OTUs numbers (Table ).

Bacterial Relative Abundance

A total of seven phyla, 14 classes, 26 orders, 43 families, and 51 bacterial genera were identified in roots, phloem, and bark of Arizona and Durango pines samples. The taxa number from phylum to family were not the same among tissues and between pine species (Table ), their frequencies varied slightly between replicates of the same tissue, among tissues of the same pine species, and between pines species. In particular, Proteobacteria was the most abundant phylum in all samples, followed by Acidobacteria, Actinobacteria, Firmicutes, Bacteroidetes, Thermi, Tenericutes, and several taxa unassigned to a phylum (Supplementary Figure S1A); among classes, Gammaproteobacteria was the most dominant, followed by Alphaproteobacteria, Acidobacteriia, Betaproteobacteria, Actinobacteria, and other in low frequencies (Figure ); and finally, among families, Enterobacteriaceae together with Acetobacteraceae and Acidobacteriaceae represented ∼80% total reads (Supplementary Figure S1B).

Distribution of relative abundance of endophytic bacterial associated with root, phloem, and bark of Durango and Arizona pines at the class level (RootPD1 = Root P. durangensis1, RootPD2 = Root P. durangensis2, PhloemPD1 = Phloem P. durangensis1, PhloemPD2 = Phloem P. durangensis2, BarkPD1 = Bark P. durangensis1, BarkPD2 = Bark P. durangensis2; RootPA1 = Root P. arizonica1, RootPA2 = Root P. arizonica2, PhloemPA1 = Phloem P. arizonica1, PhloemPA2 = Phloem P. arizonica2, BarkPA1 = Bark P. arizonica1).

A total 51 bacterial genera were recovered combining manual and QIIME assignment, among which Acetobacter (14.20%), Burkholderia (13.55%), Caulobacter (11.15%), Pseudomonas (10.27%), Ralstonia (10.03%), Bradyrhizobium (7.00%), and Methylocapsa (6.67%) were the best represented genera, followed by Rhizobium (3.64%), Providencia (3.33%), Halomonas (3.03%), Stenotrophomonas (1.76%), Serratia (1.74%), Propionibacterium (1.36%), Enterobacter (1.29%), and Mesorhizobium (1.10%); the 36 remaining genera exhibited frequencies < 1.00% of total reads (Figure and Supplementary Table S1).

Heatmap of relative abundance of endophytic bacterial genera in Durango and Arizona pines. Bacterial community distribution of the 51 bacterial genera among root, phloem, and bark samples (cophenetic correlation = 0.90) (RootPD1 = Root P. durangensis1, RootPD2 = Root P. durangensis2, PhloemPD1 = Phloem P. durangensis1, PhloemPD2 = Phloem P. durangensis2, BarkPD1 = Bark P. durangensis1, BarkPD2 = Bark P. durangensis2; RootPA1 = Root P. arizonica1, RootPA2 = Root P. arizonica2, PhloemPA1 = Phloem P. arizonica1, PhloemPA2 = Phloem P. arizonica2, BarkPA1 = Bark P. arizonica1).

α- and β-Diversity Analysis

Overall, the Good’s coverage suggested that the sampling effort for bark samples of Durango and Arizona pines was acceptable (>63%), and for phloem and root samples it was appropriate (>93%), indicating that the current sampling effort was adequate to obtain the most abundant OTUs (Table ). Rarefaction curves exhibited similar results, and these tended toward saturation in the phloem and root samples, but not in bark samples, where more sampling effort is necessary (Supplementary Figure S2). No significant differences were observed in bacterial richness and diversity among tissue samples of both pine species with Chao1 (Welch F-test: F = 4.59, p = 0.169) and Shannon index (Welch F-test: F = 10.62, p = 0.067), respectively. In contrast, the phylogenetic diversity with PD index (Welch F-test: F = 29.41, p = 0.026) and the dominance estimated with the Simpson index (Welch F-test: F = 26.75, p = 0.020) were statistically significant, suggesting that bark samples of both pine species were different to roots and phloem samples (Table ).

The bacterial communities were very similar between pine species, sharing 39 of the 51 bacterial genera identified (Figure ). Among root, phloem, and bark samples of Durango pine, 15 bacterial genera were shared, root and bark samples shared 14 genera, and four genera were present in bark and phloem samples; 11 and one bacterial genera were exclusive of bark and root, respectively, in this pine species (Figure ). In the case of Arizona pine, 12 genera were shared among roots, phloem, and bark samples, nine genera were shared exclusively between root and bark samples, three were shared between bark and phloem samples, and only one was shared between root and phloem samples; seven were found to be unique to root, one to phloem, and 11 to bark (Figure ). Finally, the comparison between pine species and gut bacteria of the bark beetle D. rhizophagus demonstrated that 13 genera were shared by these pines and the insect, 10 were exclusive of D. rhizophagus, and 38 exclusive to both pine species (Figure ; detailed information on bacterial genera available in Supplementary Table S2). The three first coordinates using the Bray–Curtis dissimilarity (Figure ), weighted UniFrac (Figure ), and unweighted UniFrac (Figure ) estimators explained 95.4, 99.68, and 56.57% of the total observed variation in each analysis, respectively. The PCoA showed significant differences between tissue samples (Bray–Curtis ANOSIM, p = 0.03; weighted UniFrac ANOSIM, p = 0.03; and unweighted UniFrac ANOSIM, p = 0.02), revealing that bark bacterial communities of both pines species were very similar, but different markedly with respect to those in phloem and root (Figures and Table ). These two last bacterial communities were similar in relative abundance, though varied in terms of richness. In the case of bacterial communities of Arizona and Durango pines and those previously reported of D. rhizophagus gut, the first three coordinates explained 93.93% of the total variation (Figure ). The PCoA showed an evident spatial segregation of communities of pine tissues and D. rhizophagus gut, with differences that were statistically significant (Bray–Curtis ANOSIM, p = 0.001). Phylogenetic inference analyses indicated there to be a high degree of taxonomic relatedness between the most abundant bacteria in the D. rhizophagus gut and their equivalent in endophytic communities (Supplementary Figure S3), that is all sequences constituted monophyletic groups at the genus level. In some cases, the topologies suggested a high taxonomic closeness at the species level (e.g., Rahnella and Pseudomonas), but others (e.g., Acinetobacter, Propionibacterium, Serratia, and Stenotrophomonas) suggested the presence of different species.

Venn diagrams showing shared and unique genera. (A) P. durangensis vs. P. arizonica, (B) D. rhizophagus vs. P. arizonica and P. durangensis, (C) P. durangensis (root vs. phloem vs. bark), and (D) P. arizonica (root vs. phloem vs. bark).

Principal coordinate analyses (PCoA) of the endophytic bacterial communities associated with root, phloem, and bark of Durango and Arizona pines samples. (A) Bray–Curtis dissimilarity, (B) weighted and, (C) unweighted UniFrac distances; and PCoA of the endophytic bacterial communities associated with root, phloem, and bark of Durango and Arizona pines and D. rhizophagus samples (D) Bray–Curtis dissimilarity. The percentages of variation explained by each axis are shown in parentheses. Orange diamonds represent the root samples. Yellow diamonds represent the phloem samples. Blue diamonds represent the bark samples. Red diamonds represent the D. rhizophagus samples (RootPD1 = Root P. durangensis1, RootPD2 = Root P. durangensis2, PhloemPD1 = Phloem P. durangensis1, PhloemPD2 = Phloem P. durangensis2, BarkPD1 = Bark P. durangensis1, BarkPD2 = Bark P. durangensis2; RootPA1 = Root P. arizonica1, RootPA2 = Root P. arizonica2, PhloemPA1 = Phloem P. arizonica1, PhloemPA2 = Phloem P. arizonica2, BarkPA1 = Bark P. arizonica1).

Functional Profiling Prediction

The PICRUSt analysis yielded average weighted NSTI values of 0.24 ± 0.07sd and 0.20 ± 0.09sd for Arizona and Durango pines samples, indicating that functional gene predictions were accurate in all samples. Among the predicted KEGG pathways at level 2, the most significant were those related to membrane transport (14.56%), energy metabolism (11.59%), carbohydrate metabolism (10.66%), amino acid metabolism (10.18%), and metabolism of cofactors and vitamins (6.68%) in terms of total pathways (Figure ). The metabolic pathways at level 3, which considered metabolic capabilities more specific, were: oxidative phosphorylation and methane metabolism (energy metabolism); glycolysis-gluconeogenesis, and metabolism of different substrates, such as starch, sucrose, pyruvate, amino sugar, and nucleotide sugars (carbohydrate metabolism); terpenoid backbone biosynthesis and prenyltransferases (metabolism of terpenoids and polyketides); glutathione metabolism, ubiquinone and other terpenoid-quinone biosynthesis, pantothenate and CoA biosynthesis (metabolism of cofactors and vitamins) (Supplementary Figure S4). The Spearman correlation of metabolic profiles showed there to be a strong positive association between endophytic and gut bacterial communities (r = 0.94; p < 0.05; Supplementary Table S3).

Heatmap depicting the PICRUSt-inferred gene relative abundance in the predicted endophytic bacterial communities of root, phloem, and bark of Arizona and Durango pines samples (cophenetic correlation = 0.98). Warm colors represent high abundances and clear colors represent low abundances (RootPD1 = Root P. durangensis1, RootPD2 = Root P. durangensis2, PhloemPD1 = Phloem P. durangensis1, PhloemPD2 = Phloem P. durangensis2, BarkPD1 = Bark P. durangensis1, BarkPD2 = Bark P. durangensis2; RootPA1 = Root P. arizonica1, RootPA2 = Root P. arizonica2, PhloemPA1 = Phloem P. arizonica1, PhloemPA2 = Phloem P. arizonica2, BarkPA1 = Bark P. arizonica1).

Discussion

In the present study, the endophytic bacterial community in roots, phloem, and bark of healthy saplings of Arizona and Durango pines was analyzed with pyrosequencing 454. Overall, Enterobacteriaceae, Acetobacteraceae, and Acidobacteriaceae were the endophytic dominant families. Previous studies using NGS technologies have demonstrated the dominance of these families (and their corresponding phyla) in needles of other mature coniferous trees, e.g., Picea engelmannii, Pinus flexilis, P. contorta, and P. radiata (Carrell and Frank, 2014; Carrell et al., 2016; Rúa et al., 2016). This information suggests that dominance of these taxa in conifers is apparently independent of the host species and age, type of tissue, geographic location, or sample processing (i.e., surface sterilization, DNA extraction protocol). Several constrains or selective factors may be responsible for this asymmetry in favor of these families, including the tree inner-environment, local interactions among different bacterial groups, bacterial metabolic capacities, and the historical association between these phyla and conifers (Hardoim et al., 2015).

At the genus level, a total of 51 bacterial genera were identified in both pine species, from which 39 genera were shared, six were exclusive to Arizona pine, and six exclusive to Durango pine. Previous studies using both culture-dependent and culture-independent techniques (e.g., DGGE, molecular cloning) have reported from one to eight bacterial genera, among which are Bacillus, Brevibacillus, Brevundimonas, Burkholderia, Cellulomonas, Kocuria, Methylobacterium, Paenibacillus, Pseudomonas, and Rahnella (Shishido et al., 1995; Pirttilä et al., 2000; Cankar et al., 2005; Izumi et al., 2008; Bal et al., 2012), representing <10% of the bacterial genera recovered in this study using pyrosequencing. These differences reflect, as it is known, the limited statistical coverage of conventional methods; however, traditional culture methods remain necessary for obtaining more detailed information about functional capacities.

The comparison between both pine species at the genus level exhibits that bacterial communities are similar in presence and abundance of taxa. It is difficult to establish a comparison at this taxonomic category with previous studies using NGS technologies, because most of them have not reported this information (Carrell et al., 2016) or there is only partial information available (Carrell and Frank, 2014, 2015; Rúa et al., 2016). However, we assume that scarce differences observed may be because of several factors as host–plant genotype, plant organs, vegetative stage, soil type, and environmental stress as was described for agricultural crops where structuring of the endophytic bacteria communities has been evident (Sturz et al., 1998; Reiter et al., 2002; Sessitsch et al., 2002; Conn and Franco, 2004).

Most of the dominant bacterial genera identified in this study (e.g., Bradyrhizobium, Burkholderia, Pseudomonas, Ralstonia, and Rhizobium) were present in all plant tissues of both pines species, whereas certain genera with frequencies <1.0% were exclusive to a particular tissue. This result indicates the lack of structuring of the dominant members of the community at genera level with respect to a specific tissue, which is the opposite of what was reported in other works (e.g., Sun et al., 2008; Reinhold-Hurek and Hurek, 2011; Garcias-Bonet et al., 2012; Xia et al., 2015; Liu et al., 2016). Whereas the tissue specificity may be explained by differences with regards to range of physical and chemical environmental conditions, such as light, oxygen concentration gradient, and toxic metabolites present in each tissue (Garcias-Bonet et al., 2012), the non-structuring may only be possible assuming that these factors in the inner-environment of conifers are similar or homogeneous in different tissues, which, although it has not been evaluated, seems to be highly unlikely. Another possible explication to be considered is that the dominant members have metabolic capacities that allow them to face and tolerate heterogeneous inner-environments and restrain community members that do not have these capacities in specific tissues. Beyond these assumptions, the structuring of these communities could be observed at finer taxonomic levels (i.e., species or strain) than those resolved by the V1–V3 16S regions in this study.

The α diversity of bacterial communities in both pine species was higher in bark than in other tissues. This higher diversity could be understood via several factors, including horizontal transmission by airborne contamination (Whipps et al., 2008), raindrops (Morris, 2002), and dissemination between plants (Vandenkoornhuyse et al., 2015). In addition, the bark is rich in organic nutrients, resulting in it being a target for many different organisms (e.g., insects, vertebrates, fungi, bacteria) (Franceschi et al., 2005), particularly in the early years of cambium activity. Regarding the lower diversity observed in roots, this might be because most coniferous species grow in suboptimal soils and climates (Moyes et al., 2016); in fact, higher diversity values have been determined in roots versus other tissues in crop plants growing in soils with greater quality and nutritional content (Albareda et al., 2006; Zarraonaindia et al., 2015). In the case of phloem in pines, diversity values obtained using culture-dependent techniques have been very low (<10 colony-forming units) (Mason et al., 2015), and to our knowledge, there are no reports of using culture-independent methods.

The PCoA analyses employing all metrics clearly separate bacterial communities of bark samples from those in root and phloem in both pine species. These differences are derived by some of the most abundant genera (Acetobacter and Methylocapsa) in bark and bacteria present or absent (Methylobacterium, Kocuria, Arthrobacter, Staphylococcus, Paracoccus, Comamonas, Prevotella, Salmonella, Shigella, Sphingomonas, Rheinheimera, Achromobacter, Capnocytophaga, Leucobacter, and Friedmaniella) in the different communities and with low abundance (<1.0% reads) (p < 0.05). Unfortunately, it is difficult to determine whether these differences are consistent among other pine or coniferous species because, until this work, there have been no other studies evaluating endospheric tissues in roots, phloem, and bark using NGS technologies. Additionally, owing to the low number of biological replicates in this study, caution in interpreting inferences about the presence/absence of particular endophytic genera in the tissues of pines is prudent.

The comparison between communities of pines and D. rhizophagus gut shows statistically significant differences (p < 0.05); however, several of the endophytic bacterial genera identified in this study apparently have taxonomic relatedness with the most abundant members of the core bacteriome of this bark beetle (Briones-Roblero et al., 2017a). While the communities of Arizona and Durango pines are more diverse in a 2:1 ratio than those previously found in D. rhizophagus, 13 bacterial genera of this insect are shared with hosts. An aspect that should be highlighted is that among shared genera, the most abundant in pines are not necessarily the most abundant in the insect gut, and vice versa, which might be a consequence of the different environmental conditions where bacteria develop (endosphere and gut). Several of these shared genera (e.g., Pseudomonas, Serratia, Rahnella) have also been found in subcortical galleries and surrounding uninfested tissues of trees colonized by other Dendroctonus bark beetles (Durand et al., 2015; Mason et al., 2015). Likewise, other shared taxa (e.g., Acinetobacter, Burkholderia, Enterobacter, Kocuria, Methylobacterium, Pantoea, Prevotella, Propionibacterium, Providencia, Pseudomonas, and Stenotrophomonas) have been reported as endophytes of pines and in numerous tissues of plants (Pirttilä et al., 2000; Strzelczyk and Li, 2000; Bal et al., 2012; Hardoim et al., 2015; Kõiv et al., 2015; Moyes et al., 2016; Yang et al., 2017).

Although it is unknown how this bark beetle acquires their symbiotic bacteria, the degree of taxonomic relatedness of bacterial genera between both assemblages suggests that at least some gut bacteria might be associated with those from pine tissues. In fact, phylogenies of the dominant bacterial genera in the gut, their equivalent endophytes, and other reference sequences of these same bacterial genera of the GenBank, formed robust monophyletic groups (boostrap value > 50%). Similarly, Hernández-García et al. (2017) found a lack of phylogenetic consistence between Dendroctonus bark beetles and their bacterial assemblages, suggesting that bacteria might be environmentally acquired. While the horizontal transmission of endophytic communities has been demonstrated in other hemimetabolous insects that feed on phloem (Lòpez-Fernàndez et al., 2017), this has not been tested in holometabolous insects as bark beetles. Specific experimental studies under controlled laboratory conditions are necessary to determine this apparent horizontal transmission in bark beetles.

Lastly, the positive association found between the metabolism of endophytic and bark beetle gut bacterial communities suggests that certain predicted functional pathways are shared by both communities (e.g., amino acids, carbohydrates, biosynthesis of secondary metabolites, energy and cofactors, vitamins metabolism, xenobiotic biodegradation). Some of these metabolic routes have also been reported in bacterial communities associated with other phloem- and wood-feeding insects as well as other Dendroctonus species (Supplementary Figure S5) using PICRUSt or shotgun metagenomics (Scully et al., 2013; Berasategui et al., 2016; Hernandez-García et al., unpublished data). In addition, certain particular metabolic functions (e.g., hydrolysis of lipids, starch, esters, xylan, and cellulose, as well as xenobiotic biodegradation) have been demonstrated in bacteria isolated from the gut of particular Dendroctonus bark beetles (Morales-Jiménez et al., 2012, 2013; Hu et al., 2014; Cano-Ramírez et al., 2016; Briones-Roblero et al., 2017b). Likewise, several of these metabolic pathways (oxidative phosphorylation, xenobiotic degradation, methane, glutathione, amino sugar and nucleotide sugar metabolism, and carbohydrates metabolism) have been described in endophytic bacterial communities of many plants using PICRUSt, and in some cases confirmed with bacterial isolates (Van Aken et al., 2004; Parmentier et al., 2011; Sheibani-Tezerji et al., 2015; Ruiz-Pérez et al., 2016; Sánchez-López et al., 2017).

In brief, this is the first study that investigates the diversity and structure of the endophytic bacterial community in roots, phloem, and bark of Arizona and Durango pines. Our findings show there to be a bacterial community similar to those reported for other pine species and coniferous genera. Our results also demonstrate that few dominant members of these communities are shared among different types of tissues and some low-abundance taxa are exclusive to particular tissues. Several genera of this endophytic community showed there to be taxonomic relatedness with gut-dominant members of the bark beetle, D. rhizophagus, suggesting that some could be the same species. However, further evidence is necessary to clarify if they are the same species or not, and corroborate their mode of transmission. Lastly, our results indicate there is a strong metabolic association at the community level, although large differences in ecological and functional traits could be observed at the species level within the same genera or families, and as a consequence, they can display different ecological responses depending on the habitat or environment. Future research based on “omics” technologies (e.g., genomics, transcriptomics, proteomics, metabolomics), as well as the assessment of metabolic capacities in microorganisms isolated from these pines will be needed to reach a better understanding of the pine-bacteria relationship.

Author Contributions

Conception of the work: RG-E, FR-O, and GZ; funding acquisition: FR-O and GZ; analysis and interpretation of data for the work: RG-E, CB-R, RP-M, FR-O, and GZ; drafting the work: RG-E, FR-O, and GZ; revising it critically for important intellectual content: RG-E, CB-R, FR-O, and GZ. All authors approved the manuscript.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Table 1

Sample ID and taxa number identified from phylum to genus.

Host	Tissue/biological replicate	Sample ID	Taxa number at level:
			Phylum	Class	Order	Family	Genus
Durango pine	root1	RootPD1	4	7	9	15	16
	root2	RootPD2	7	11	16	23	26
	phloem1	PhloemPD1	5	9	14	20	18
	phloem2	PhloemPD2	7	9	10	14	11
	bark1	BarkPD1	6	13	24	38	42
	bark2	BarkPD2	4	7	12	17	12
Arizona pine	root1	RootPA1	6	12	17	25	25
	root2	RootPA2	5	10	16	25	21
	phloem1	PhloemPA1	6	9	13	20	16
	phloem2	PhloemPA2	2	5	9	13	10
	bark1	BarkPA1	6	11	20	32	36
	bark2	BarkPA2	–	–	–	–	–

Table 2

Richness and α diversity indices for root, phloem, and bark samples of Durango and Arizona pines.

Host	Tissue	Good’s coverage (%)a	Observed OTUs97%,a	Chao1a	Simpsona	Shannona	PDa,b
Durango pine	Root	94	25	116	0.56	2.1	1.77
	Phloem	96	20	46	0.48	1.79	1.21
	Bark	63	144	529	0.96	6.25	9.38
Arizona pine	Root	93	29	111	0.58	2.25	2.56
	Phloem	94.5	21	104	0.48	1.77	1.11
	Bark	74	97	443	0.91	4.89	7.17