diversity Thrips Microbiome Study in Commercial Avocado ( Persea americana Mill.) from Northwest Colombian Andes (Antioquia, Colombia) Shows the Presence of Wolbachia , Ehrlichia , Enterobacter

: Microbiota associated with insects play several important roles in their host, including protection against pathogens, provision of nutrition, and survival in hostile environments. The aim of this work was to identify the bacterial community found in avocado thrips from Northwestern Colombia (Antioquia department) in order to ﬁnd isolates for potential biocontrol purposes. Culture-dependent methods based on 16S rRNA and gyrase B gene sequencing in 42 bacterial isolates allowed the identiﬁcation of the genera Bacillus , Serratia , Moraxella , Pantoea , and Sphingomonas . Microbial diversity detected with the temperature gradient gel electrophoresis (TGGE) technique on three morphotypes of thrips, named brown ( Scirtothrips hansoni) , black ( Frankliniella panamensis ), and pale ( Frankliniella sp.), showed a low bacterial community density (Shannon–Wiener index = 1480, p > 0.05) with signiﬁcant differences among morphotypes (R = 0.7877, p = 0.0004). Results obtained with Illumina sequencing on the V1–V2 hypervariable region of the subunit 16S rRNA showed a predominant sequence in the brown morphotype ( Scirtothrips hansoni ) that belongs to the genus Wolbachia . The 16S amplicon analyses were extended to more samples and higher resolution using the V4–V5 hypervariable region. The results showed six additional bacteria phyla, conﬁrming the previous observation for the dominant bacterial groups made in S. hansoni and the detection of the alternation of highly predominant genera among these thrips. Our results demonstrate that endosymbiont such as Wolbachia sp. are part of the microbiota of these pests, thereby indicating the possibility of employing this type of bacterium to improve the management of avocado thrips globally.


Introduction
Avocado (Persea americana Mill.) is an evergreen, fruit-yielding tree. The production of this fruit (about six million tons) has significantly increased, and approximately 62% of the global crop is harvested in five main countries, namely, Mexico (33%), the Dominican Republic (10.5%), Peru (7.8%), Indonesia (5.7%) and Colombia (5.1%) [1] However, this crop  Figure 1). The collections were carried out with oral aspirators [28] during the months of March and August 2017. The insects were transported live to the Cellular and Molecular Biology laboratory at Universidad Nacional de Colombia, headquartered in Medellin, and stored at 4 • C for later processing. A total of 731 thrips were collected, and 381 thrips were used for culture-dependent analysis and 350 thrips were used for culture-independent analysis.

Isolation and Purification of Bacteria from Thrips
Aerobic and cultivable bacteria were isolated from 381 adult specimens that were immediately processed after collection (Table S1). The samples were washed with 0.5% hypochlorite for 10 s then immersed in 70% ethanol for 10 s and washed three times in sterile water for 10 s to remove excess dust and exogenous bacteria. Each sample included 1 or 25 thrips that were homogenized to produce a pool used for serial dilutions in a 1 X Phosphate Saline Buffer (PBS) and further plated for bacterial growth on Luria-Bertani (LB) agar (BD Difco, East Rutherford, NJ, USA) and medium R2A agar (Merck KGaA, Darmstadt, Germany; [29,30]). Petri dishes were aerobically incubated at 25 • C for 24 to 48 h for LB and 7 days for R2A. LB and R2A media were chosen as non-selective media to promote the growth of a greater diversity of cultivable microorganisms. The total number of colonyforming units (CFUs) was estimated for each homogenate, and different colonies were selected based on the macroscopic (colony size, pigmentation, elevation) and microscopic characteristics by Gram stain. These were isolated from the different media plates and sub-cultured to obtain pure bacterial cultures. All isolates were cryopreserved in 20% glycerol at -80 • C. Comparisons between the mean number of CFUs between the Hass and Fuerte varieties were made using a t-test using the Past 3.0 software [31]. The t-test was estimated after performing a Kolmogorov-Smirnov normality test and a variance homogeneity test.

Extraction of Genomic DNA and PCR of the Transcribed Internal Spacer-ITS, 16S rRNA and gyrB Genes
For the amplification of ribosomal genes, bacterial genomic DNA was extracted from the selected bacteria using a heat lysis protocol [32] from freshly cultured cells. This DNA was later used as a template for PCR assays and was diluted at 1:10 in sterile water. Template DNA from all isolates was initially evaluated by amplifying the inter-ribosomal region (ITS) using primers L1 (5 -CAAGGCATCCACCGT-3 ) and G1 (5 -GAAGTCGTAACAAGG-3 ) [33], as previously described by Moreno et al. [34]. The amplified bands were run in polyacrylamide gels (8%) and then the software GelCompar II (Applied Biosystems Maths, Belgium) was used to construct a dendrogram, based on the obtained band patterns, using Pearson correlation distance matrices and complete-linkage grouping methods [35]. Fifty percent or more similarity values between ITS standards were established as criteria for subsequent molecular identification.
The DNA isolated from selected colonies was used to amplify the 16S rRNA gene (1.5 kbp), using the primers Eubac 27 F (5 -AGAGTTTGATCCTGGCTCAG-3 ) and 1492R (5 -GGTTACCTTGTTACGACTT-3 ) according to Moreno et al. [34]. This gene was further sequenced for the identification of each isolate affiliation. Additionally, bacterial DNA was extracted for amplification of the gyrB gene using the UltraClean Microbial DNA Isolation Kit (MOBIO) following the manufacturer's instructions. Primers UP1 (5 -AGC AGG GTACGGATGTGCGAGCCRTCNACRTCNGCRTCNGTCAT-3 ) and UP2R (5 -GAAGTC ATCATGACCGTTCTGCAYGCNGGNGGNAARTTYGA-3 ) [36] were used in a PCR reaction to amplify a conserved region of approximately 1.26 kbp of the gyrB gene [37]. This locus was also used to corroborate the identification of each affiliation of the isolates. For this gene, the following thermal profile was used: 1 cycle of 94 • C for 5 min; 30 cycles of: 94 • C for 1 min, 53 • C for 1 min, 72 • C for 90 s; and a final cycle of 72 • C for 7 min. DNA positive (DNA from a pure culture of Bacillus cereus) and negative (ultrapure water) controls were routinely included in all PCR reactions. Amplification products were verified by visualization in 1% agarose gel and purified using Illustra GFX PCR DNA and a Gel Band Purification Kit (GE Healthcare, UK), and the double-stranded DNA was sequenced in both directions using the ABI PRISM 3700 DNA analyzer service of Applied Biosystems.

Temperature Gradient Gel Electrophoresis (TGGE)
The TGGE protocol was used to monitor the dynamic changes in the microbe population in the three morphotypes. Total DNA samples were quantified using a ND-100 nanodrop Thermo Scientific spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA). Subsequently, these samples were subjected to amplification of the 16S rRNA gene (variable regions V3-V6). These samples were amplified with forward primer 341 F (5 -CCTGCAGGAGGCAGCAG-3 ), with an additional GC termination at the 5 end and the reverse primer 907R (5 -CCCTGACGTGTTATTCAATTCY-3 ). PCR reaction conditions were implemented as previously described for 16S rRNA gene amplification [39].
Gels were stained with 1X GelGreen ® Nucleic Acid Gel Stain for 30 min and the obtained image was digitized using Gel-Compare II 6.6 Software (http://www.appliedmaths.com accessed on 8 April 2019). A cluster analysis was constructed through Pearson correlation and complete-linkage clustering [35], to detect differences in diversity and abundance profiles amongst samples. Bands of interest were cut, and the DNA was eluted in 30 µL of ultrapure water. Five microliters of the DNA elution were re-amplified using the same primers, 341 F (without GC tail) and 907R. Additionally, some bands were selected to be cloned with the CloneJET ® PCR Cloning Kit (Thermo Fisher Scientific, Waltham, MA, USA) according to the manufacturer's instructions. Bands and clones that were successfully re-amplified were sent for sequencing using the methods and procedures described above.

Sequencing and Phylogenetic Analysis
Sequences from PCR products of the 16S rRNA gene, gyrB gene of bacterial isolates and TGGE eluted bands were edited and cured using CLC Genomics Workbench software (Version 7.0.3) (CLC Bio, Aarhus, Denmark). Consensus sequences were contrasted with the GenBank database using the BLASTn algorithm. 16S rRNA sequences were also compared with Ribosomal Database Project (RDP) [40,41]. The homologous sequences from the Genbank database were aligned together with sequences obtained here using the MUSCLE function [42] found in MEGA X [43]. The sequences obtained here were further sent to the GenBank database.
The phylogenetic affiliations amongst sequences 16S rRNA and gyrase were analyzed using a Bayesian tree constructed in MrBayes (Version 3.2.0) [44] with an MCMC chain using 10 million generations and four chains for every 1,000,000 generations. Verification and validation of recombination events and the presence of chimeras were carried out with RDP2 software [41]. Methanococcus aeolicus (DQ136171) and Bacteroides caccae (NZ_CP022412) were used as outgroups for phylogenetic analysis of the 16S rRNA and gyrB genes, respectively. For the 16S rDNA partial sequences obtained by TGGE, a dendrogram was constructed with the neighbor-joining (NJ) and Kimura 2-parameter methods as evolutionary distances; finally, bootstrap analysis was undertaken using 5000 iterations with MEGA X [43].

Analysis of the Bacterial Community by PCR-TGGE
A binary matrix was constructed from the presence/absence of bands evidenced in the TGGE gel in the program PAST (Version 3.23). Principal component analysis (PCA) was used to determine differences amongst samples, using 10000 bootstraps and confidence intervals for the execution that represent 95% of the data [30]. Moreover, a grouping was performed with the unweighted pair group algorithm (UPGMA) and the Dice similarity index, and with 10,000 bootstraps. With the same binary matrix, a unidirectional analysis of similarities (ANOSIM) was performed using the Dice index to estimate the similarity of values between bacterial morphotypes of thrips of the bands obtained by TGGE [45].

16S Amplicon Next Generation Sequencing (NGS)
The quality of the DNA of each group of thrips separated by morphotype (n = 10), obtained as previously described, was verified by spectrophotometry in a Nanodrop ® 2000 (Thermo Fisher Scientific, USA). PCR amplification was performed using primers 27F (5 -AGAGTTTGATCMTGGCTCA-3 ) and 338R (5 -GCTGCCTCCCGTAGGAGT-3 ) [46] targeting hypervariable region V1 and V2 of the 16S rRNA gene, initially used to analyze three samples belonging to S. hansoni (brown morphotype) coming from two farms: LV (La Vega) and VFm. Further on, as is standard in studies of environmental microbiomes of the Earth Microbiome Project (EMP https://earthmicrobiome.org/protocols-and-standards/16s/ accessed on 17 April 2019), subsequently, the primers 515F (5 -GTGYCAGCMGCCGCGGTAA-3 ) [47] and 806R (5 -GGACTACNVGGGTWTCTAAT-3 ) [48] targeting the hypervariable region V4-V5 of 16S rRNA genes were used in the following samples as they could provide additional resolution inside the taxonomic groups found in the initial samples and could also detect Archaea. The samples were: three samples of S. hansoni (brown morphotype) were analyzed and taken from samples of farms: LV2 (La Vega), AASS1 (San Sebastián 1 = sampling period 1) and AASS2 (San Sebastián, sampling period 2); three Frankliniella samples (pale morphotype) from the farms: BH1 (Barahonda), EP1 (El Paraiso), and LS1 (Luisines); and a sample of F. panamensis (dark morphotype) from the farm LVn (La Vega). The PCR products were used for library construction using Illumina MiSeq platform 2 × 250 bp. The obtained sequences were paired and processed following the QIIME 1.9.1 pipeline for the V1-V2 region and QIIME 2 for the V4-V5 region [49]. Chimeras were removed and raw sequences were quality filtered, clustered into operational taxonomic units (OTUs) and taxonomically classified using the "Naive Bayesian Classifier" of the Ribosomal Database Project (RDP) classifier against the SILVA 132 database [50].
Sample comparisons were performed in terms of relative abundance, and alpha and beta diversity, using R scripts developed by Pérez-Jaramillo et al. [51] following the procedure of Castañeda-Monsalve et al. [27]. Results were plotted using the Phyloseq R package (V 1.24.2) [52]. Comparisons of bacterial diversities based on the bacterial genera dataset were performed with a chi-square test performed with R and significant p-values were considered after a Bonferroni correction (α = 0.05/28 = 0.0018). Observed OTUs, Shannon, Simpson, and Chao1 diversity metrics were also obtained. One-way ANOVA and Tukey HSD were used to test for significance. Additionally, a Venn diagram was calculated using the DISPLAYR web platform (https://www.displayr.com/ accessed on 14 January 2022) to determine the number of specific and shared genera among thrips morphotypes. These analyses were performed on the set of sequences of hypervariable region V4-V5 16S rRNA gene amplicons. For beta diversity calculations, OTU table normalizations were undertaken as described by Pérez-Jaramillo et al. [51]. Differences were visualized using principal coordinate analysis (PCoA), and constrained analysis of principal coordinates (CAP) [53]. Permutational multivariate analysis of variance (PERMANOVA) was used for hypothesis testing; differences between factors were considered to be significant if the likelihood of the observed statistics was <5% (p < 0.005).

Isolation and Culture of Microbial Isolates
Bacterial growth was followed in 381 adult thrips that were processed in groups of 15 insect pools containing from 1 to 25 specimens (Table S1). Mean bacterial count ranged from 2.6 × 10 4 to 3.64 × 10 4 colony forming units (CFU/mL) for Hass and Fuerte varieties, respectively. No significant differences were found in these counts for the two avocado varieties studied (Hass vs. Fuerte) (Table S1).
After examining the macro and microscopic characteristics with Gram staining, 42 isolates were selected and further used for molecular characterizations based on an analysis pattern of the transcribed internal spacer region (ITS) ( Figure S1). Profiles of ITS region showed significant similarities between isolates, with common bands between 700-900 bp and 400-600 bp. The dendrogram differentiated 29 clusters with a similarity percentage of 50%. Representative isolates of each cluster were further selected for 16S rRNA gene sequencing. No differences were found in bacterial isolates between Hass and Fuerte varieties, meaning that both share similar bacterial isolates (Table S1, Figures S1 and S2).

Phylogenetic Analysis with Sequences of the 16S rRNA and gyrB Genes
Sequence-based molecular identification of the 16S rRNA and gyrB genes showed a high level of correspondence with sequences deposited in the NCBI at the species level and between different taxonomic divisions (Table S2). The Bayesian tree obtained with 16S rDNA sequences ( Figure 2a) and gyrB (Figure 2b) produced two clusters; the first corresponded to the Bacillus genus and the second to Sphingomonas, Moraxella, Pantoea, and Serratia genera. Similarly, phylogenetic results showed that the Erwinia aphidicola (BFo1) strain was found as a symbiont in F. occidentalis (Chanbusarakum and Ullman, [14]; Facey et al. [18]) and this species has 97% identity with our sample T8H6 ( Figure 3a). Furthermore, the sample T3H1_1 exhibits a 99% identity percentage with Pantoea agglomerans also reported in other thrips populations and thrips species such as F. occidentalis (De Vries et al. [54]) and T. tabaci (De Vries et al. [19]).

Phylogenetic Analysis with Sequences of the 16S rRNA and gyrB Genes
Sequence-based molecular identification of the 16S rRNA and gyrB genes showed a high level of correspondence with sequences deposited in the NCBI at the species level and between different taxonomic divisions (Table S2). The Bayesian tree obtained with 16S rDNA sequences ( Figure 2a) and gyrB (Figure 2b) produced two clusters; the first corresponded to the Bacillus genus and the second to Sphingomonas, Moraxella, Pantoea, and Serratia genera. Similarly, phylogenetic results showed that the Erwinia aphidicola (BFo1) strain was found as a symbiont in F. occidentalis (Chanbusarakum and Ullman, [14]; Facey et al. [18]) and this species has 97% identity with our sample T8H6 ( Figure 3a). Furthermore, the sample T3H1_1 exhibits a 99% identity percentage with Pantoea agglomerans also reported in other thrips populations and thrips species such as F. occidentalis (De Vries et al. [54]) and T. tabaci (De Vries et al. [19]).

PCR-TGGE Analysis
In total, 350 thrips were processed to obtain 14 pools of DNA that were further used for the PCR-TGGE analysis (Table S3). This technique allowed the differentiation of fragment profiles of the gene 16S rRNA of approximately 585 bp ( Figure S3), assuming that each band (with a different pattern) represented a different taxonomic unit (OTU). From this gel, 68 unique and common DNA bands were selected based on their intensity and were then sequenced and identified according to their relationship with other reported sequences from GenBank. TGGE patterns observed for each sample were similar between replicates ( Figure 3 and Figure S3).
After analyzing band patterns for each profile with GelCompar II ® software, two groups were observed with similarity rates of 50% and 60% found in the UPGMA dendrogram. The groups of profiles obtained here clustered according to the three morphotypes of thrips studied, as shown in Figure 3

PCR-TGGE Analysis
In total, 350 thrips were processed to obtain 14 pools of DNA that were further used for the PCR-TGGE analysis (Table S3). This technique allowed the differentiation of fragment profiles of the gene 16S rRNA of approximately 585 bp ( Figure S3), assuming that each band (with a different pattern) represented a different taxonomic unit (OTU). From this gel, 68 unique and common DNA bands were selected based on their intensity and were then sequenced and identified according to their relationship with other re- Based on the migration and band intensity obtained from the TGGE profile, the bacterial diversity and richness amongst the samples were estimated. The Shannon-Wiener index was calculated, producing values of 1480 and 15.6404, respectively. Both estimators were low for the thrips studied.
Results obtained with ANOSIM analysis showed significant differences between the band patterns found in the three morphotypes of thrips (R = 0.7877, p = 0.0004) (Figure 3 and Figure S3). Wolbachia and non-cultivable Enterobacteriaceae bacteria prevailed in thrips of the brown morphotype (S. hansoni), non-cultivable Firmicutes bacteria were abundant in the pale morphotype (Frankliniella), and Pantoea and Anaplasma were also abundant in dark morphotype (F. panamensis).

Phylogenetic Affiliation of the Sequences Obtained from the TGGE Bands
For phylogenetic analysis, several bands were selected from profiles representing a bacterial operating taxonomic unit (OTU) (Figure 4). Forty-eight bands and seven clones from a total of sixty-eight bands produced good quality sequences that were then assembled, aligned, and identified with high percentages of similarity (98-100%). These bands correspond to specific bacterial groups sequenced with the 16S rDNA based on public databases (Table S4). The affiliations found showed that the phylum Proteobacteria was dominant, with Wolbachia being the most dominant genus (99%) in all samples; Burkholderia, Pantoea, and Ralstonia were also found (Table S4). In the phylum Firmicutes, the Bacillus genus and non-cultivable Firmicutes were found, and, for the phylum Actinobacteria, the genus Cutibacterium was found (Table S4).
Based on the migration and band intensity obtained from the TGGE profile bacterial diversity and richness amongst the samples were estimated. The S non-Wiener index was calculated, producing values of 1480 and 15.6404, respecti Both estimators were low for the thrips studied.
Results obtained with ANOSIM analysis showed significant differences betwee band patterns found in the three morphotypes of thrips (R = 0.7877, p = 0.0004) (Fig  and Figure S3). Wolbachia and non-cultivable Enterobacteriaceae bacteria prevaile thrips of the brown morphotype (S. hansoni), non-cultivable Firmicutes bacteria abundant in the pale morphotype (Frankliniella), and Pantoea and Anaplasma were abundant in dark morphotype (F. panamensis).

Phylogenetic Affiliation of the Sequences Obtained from the TGGE Bands
For phylogenetic analysis, several bands were selected from profiles represent bacterial operating taxonomic unit (OTU) (Figure 4). Forty-eight bands and seven c from a total of sixty-eight bands produced good quality sequences that were the sembled, aligned, and identified with high percentages of similarity (98-100%). T bands correspond to specific bacterial groups sequenced with the 16S rDNA base public databases (Table S4). The affiliations found showed that the phylum Proteob ria was dominant, with Wolbachia being the most dominant genus (99%) in all sam Burkholderia, Pantoea, and Ralstonia were also found (Table S4). In the phylum Firmic the Bacillus genus and non-cultivable Firmicutes were found, and, for the phylum tinobacteria, the genus Cutibacterium was found (Table S4).  OTUs with readings representing less than 0.0001% of the total depth of the samples were filtered to eliminate so-called singletons as possible technical artifacts or low-frequency ribotypes that do not contribute significantly to the composition of the community. After filtering through rarefaction analysis, 24.773 reads remained for the V1-V2 region ( Figure S5), assigned to 363 OTUs, and, for the V4-V5 region, 13.000 reads ( Figure S6) were assigned to 627 OTUs.
With the V1-V2 region, two genera with a threshold >2% were identified in a sample of S. hansoni (brown morphotype). These genera were Wolbachia (93.98%) and Rosenbergiella (1.84%), while in the other samples of the same morphotypes Wolbachia prevailed with abundances of 95.89% and 92.85%, respectively (Figure 5c). Our results indicate a prevalence of the Wolbachia genus in samples of the S. hansoni species.
With the V4-V5 region, at the genus level, 28 genera were identified with a threshold >2%. This result suggests a great variation in the microbial communities amongst the three morphotypes ( Figure 5d, Table S5). The most abundant and common genera found were: Wolbachia (except for a pale EP1 sample), Ralstonia, Phaseolibacter, and Enterobacter (mainly in the pale morphotype), and Streptococcus and Sphingomonas (mainly in the pale and brown morphotypes). Furthermore, the Halorubrum genus (Archaea) was present in both brown and pale morphotypes with abundances of 0.02%. Additionally, a species-specific association was found between S. hansoni (brown morphotype) and Wolbachia genus, in addition to Ehrlichia in F. panamensis (dark morphotype) and Enterobacter and Rickettsia genera in Frankliniella (pale morphotype).
Statistical chi-square analysis showed that Wolbachia, Sphingomonas, Parvibaculum, and Deinococcus-Thermus genera were significantly more abundant in S. hansoni (brown morphotype); Thalassospira, Enterobacter, and other bacteria in Frankliniella (pale morphotype); and, finally, Ehrlichia genus in F. panamensis (dark morphotype) (Table S5). Additionally, in the Venn diagram, it was found that three OTUs were shared for the three morphotypes of thrips belonging to the Enterobacter, Ralstonia, and Pseudomonas genera (Figure 5e, Table S5).

Bacterial Diversity
For the region V1-V2, alpha diversity analysis was performed after rarefaction, taking all samples to a depth of 24,773 sequences/sample; rarefaction curves reached a plateau between 22,000 to 25,000 reads. For the region V4-V5, alpha diversity analysis was performed after rarefaction, taking all samples to a depth of 13,000 sequences/sample; rarefaction curves reached a plateau between 12,000 to 13,000 reads.
Shannon and Simpson indices for the V1-V2 region were estimated to assess bacterial diversity at the morphotype level in S. hansoni (brown morphotype) (Figure 6a). Low diversity was observed, and it was not possible to perform an ANOVA since the three samples evaluated belonged to the same brown morphotype (S. hansoni) (single level). The estimated diversities for S. hansoni (brown morphotype) samples from LV, VFm, and VFp with the V1-V2 region were as follows: Shannon = 0.74 +/−0.11, Simpson = 0.19 +/−0.03. These indices were estimated in two farms: La Vega (LV) located in the municipality of La Ceja and Villa Fátima (VFm and VFp) located in the municipality of Guarne (Figure 6a).
Similarly, Shannon and Simpson indices for the V4-V5 region were evaluated to deepen the assessment of the bacterial diversity amongst S. hansoni (brown morphotype) and to expand the analyses in Frankliniella (pale morphotype) and F. panamensis (dark morphotype). Statistical differences were observed amongst them as the diversity was the highest in S. hansoni (brown morphotype) and the lowest in the dark morphotype. Alpha diversity indices estimated for the V4-V5 region showed that Shannon and Simpson indices were the highest for the brown morphotype (Shannon mean = 1.79 +/−1.05, Simpson mean = 0.56 +/−0.23) and the lowest for the dark morphotype (Shannon mean = 0.27, Simpson mean = 0.08) (Figure 6b). It is important to mention that, for the dark morphotype, only one sample was considered. This may reflect the estimates of low bacterial communities obtained for this specimen. ANOVA results based on the differences in these three indices between the three morphotypes were not significant (Table S6). This result may be due to the normalization of the samples since the AASS2 sample presented 13.000 reads (sample with the lowest reads) and this sample represented the threshold used for this analysis.
Principal coordinate analysis (PCoA) obtained with Bray-Curtis distances was only performed to compare the bacterial communities identified on the three morphotypes with the V4-V5 region. This analysis showed that the variable "Morphotype" explained 49.7% of the β-diversity of the microbial communities associated with avocado trips (PERMANOVA, p = 0.05, Table S7) (Figure 6c). Principal coordinate analysis (PCoA) obtained with Bray-Curtis distances was only performed to compare the bacterial communities identified on the three morphotypes with the V4-V5 region. This analysis showed that the variable "Morphotype" explained 49.7% of the β-diversity of the microbial communities associated with avocado trips (PERMANOVA, p = 0.05, Table S7) (Figure 6c).

Discussion
In this research, we reveal the composition of the thrips microbiota using a culturedependent approach with a taxonomic identification carried out with sequencing of the 16S rRNA and gyrase genes combined with culture-independent techniques supported by a next-generation sequencing analysis. A low diversity was observed both at the culturedependent and culture-independent level and differences were observed between the three morphotypes evaluated and identified as S. hansoni (brown morphotype), F. panamensis (dark morphotype), and Frankliniella (pale morphotype).
With culture-dependent approaches, Proteobacteria and Firmicutes were identified as the predominant phyla. Although the culture-independent techniques showed the presence of the previously found phyla Proteobacteria and Firmicutes, plus Actinobacteria by TGGE, high-throughput sequencing of the 16S rDNA extended the detection of the microbial community to six phyla, mainly composed of Proteobacteria, with a minor fraction shared between representatives of Cyanobacteria, Actinobacteria, Deinococcus-Thermus, Bacteroidetes, and Firmicutes. With all methodological strategies used, Proteobacteria appeared as the dominant phylum in the thrips morphotypes. Although analyses of the microbiota of the order Thysanoptera are limited, especially in those that attack avocado, recent studies confirmed that Proteobacteria is the most abundant in the bacterial communities of other thrips species, such as S. dorsalis, Thrips palmi, and Hoplothrips carpathicus. Proteobacteria comprise~46 to~60% of the bacteria detected in thrips [7,11,25,55]. In contrast, Firmicutes have been found in low abundance and comprise no more than 10% [25].
Through culture-dependent methods we obtained 42 bacterial isolates according to their in vitro morphological characterization, which was selected for 16S rDNA sequencing.
Additionally, 25 strains from this group were used for the analysis of the gyrB gene sequence to corroborate the identity of bacteria with both genes (Table S1, Figure 2). The results revealed the presence of Bacillus, Serratia, Moraxella, Pantoea, and Sphingomonas genera (Table S2, Figure 3a,b). The Gram-negative Pantoea and Serratia genera have been previously identified in Frankliniella occidentalis (Pergande) and Thrips tabaci (Lindeman) associated with other crops, such as flowers and tobacco [14,16,19]. The main Gram-negative bacteria previously reported in these two species of thrips were the genera Agrobacterium, Erwinia, Methylobacterium, Pseudomonas, and Serratia [7,14,17,19].
In more recent studies, symbiotic relationships between thrips and Erwinia have been described [14,19]. Erwinia was detected in several countries and in different crops, such as beans, cucumbers, and chrysanthemums, suggesting that this genus is transmitted horizontally. Specifically, two different types of bacteria, BFo-1 and BFo-2, were identified in two populations of F. occidentalis [14]. Genome sequencing revealed that BFo-1 is more closely related to species of the Erwinia genus and is a close relative of Erwinia aphidicola [18]. Furthermore, the BFo-2 type bacterium is an enterobacterial microorganism and is not phylogenetically close to other currently known bacterial species. The BFo-2 strain detected here shares a common ancestor with Pantoea and Erwinia genera [18]. The results obtained in the present work showed a value of 97% of BFo-1 with our sample T8H6 (Figure 2a), whereas the sample T3H1_1 had a 99% identity percentage with Pantoea agglomerans. Moreover, bands obtained by culture-independent methods showed the presence of noncultivable bacteria sequences that are closely related to BFo-2 ( Figure 5), suggesting that avocado thrips harbor both bacteria also identified in F. occidentalis [14,18]. All these data suggest that these species share some common microbiota regardless of the species or the crop that they live on.
Pantoea is an interesting Gram-negative genus detected in the present study. This genus has been extensively analyzed in various insects [56] and can degrade herbicides [57] and produce antimicrobials [58]. Most importantly, Pantoea is known to be a symbiote of various insects, as it helps to degrade cellulose in plants (Suen et al., 2010). Other Gramnegative species that were identified in this study include Sphingomonas olei, M. osloensis, and S. liquefaciens. These bacteria are capable of degrading aliphatic hydrocarbons [59], molluscicides [60], and insecticides [61].
Regarding the Gram-positive bacteria found here, several Bacillus strains were identified. Identification of Bacillus spp. using the 16S rRNA gene did not allow adequate identification at species level due to the high percentage of gene sequence similarities amongst closely related species [62]. However, identification using the gyrase gene gave an approximation of the possible Bacillus species, suggesting the presence of B. safensis, B. cereus, B. thuringiensis, B. velezensis, and B. amyloliquefaciens. These bacteria have the potential to be biocontrollers of plant pathogens, since they can produce lipopeptides [63]. B. amyloliquefaciens, B. cereus, B. pumilus, B. velezensis, and B. thurigiensis, amongst other species, can produce a significant reduction in the incidence of diseases in various host plants [63][64][65].
Regarding the culture-independent analysis, differences in microbiota composition and diversity were further evaluated in three different thrips. This analysis involved the extraction of biological material from the entire insect, which means that allochthonous and indigenous microbiota could be included in the DNA samples. Data showed that the three morphotypes shared phylotypes, but with stark differences between them (Figures 4 and 5). This was additionally evidenced by a principal component analysis from the binary matrix of the banding profiles of the TGGE, which also showed significant differences between the three morphotypes ( Figure 4) due to results obtained with the similarity indexes and the grouping UPGMA ( Figure S3) and one-way ANOSIM. In culture-independent assays, and TGGE and metataxonomic methods, Proteobacteria were found most abundantly, particularly in S. hansoni (brown morphotype), a prevalence of the Wolbachia genus, and non-cultivable Enterobacteriaceae. Non-cultivable Firmicutes bacteria were found in Frankliniella (pale morphotype), and Pantoea and Anaplasma genera in F. panamensis (dark morphotype). Other genera, such as Burkholderia, Ralstonia, Bacillus, Cutibacterium, and Heliothrips haemorrhoidalis symbiotes were mainly present in Frankliniella (pale morphotype). These kinds of symbiotes have been previously described in other insects such as H. haemorrhoidalis, Telenomus tridentatus, and Aedes spidentatus [66][67][68].
Results obtained by high-throughput sequencing showed that Proteobacteria was the most abundant phylum found in the three morphotypes with a relative abundance >95 for the V1-V2 and V4-V5 regions of the 16S rRNA gene. This result corroborated previous studies of other species of thrips [7,11,25] and other insects [69,70], demonstrating the prevalence of this phylum in insects. Similarly, the bacterial communities detected with the V4-V5 region showed the presence of Firmicutes, Actinobacteria, Cyanobacteria, Bacteroidetes, and Deinococcus-Thermus: these phyla have also been previously reported in other species of thrips [7,11,25,71].
Furthermore, the relative abundance of the different genera detected here showed the prevalence of Wolbachia in S. hansoni (brown morphotype) of >92% with the V1-V2 region and >65% with the V4-V5 region, whereas the prevalence of this endosymbiont was lower for Frankliniella (pale morphotype), ranging from 0 to 2.43%, and F. panamensis (dark morphotype), equal to 0.16%.
Wolbachia has been previously identified in other species of thrips. For example, a study by Arakaki et al. [21] found Wolbachia supergroup B in Franklinothrips vespiformis microbiota. Moreover, Dickey et al. [7] detected this endosymbiont in S. dorsalis. Furthermore, Saurav et al. [24], observed that the female and male proportions found in the natural populations of T. palmi were 3:1 and identified Wolbachia in the abdomen of the insect. They suggested that differences in the proportions of the sexes were due to an incompatibility cytoplasmatic effect of this bacteria between genera in this insect. Another study by Kaczmarczyk et al. [25] found that, in H. carpathicus, the Wolbachia genus was present in 70% of the pupae and in 56% of adults. Ambika and Rajagopal [72] also found that Wolbachia was associated with Plicothrips apicalis in the abdomen of their larvae and adults. Additionally, observations made in this study revealed a higher proportion of females relative to males in nature, particularly for Frankliniella, whereas, for S. hansoni, no males were collected in the field. The presence of Wolbachia may be related to the differences in male and female proportions in these two species. Moreover, reproduction by parthenogenesis has been reported in thrips [21] and this may be due to an effect of Wolbachia. Detection of this bacteria is relevant in this study as further investigations should be undertaken to determine its prevalence in avocado thrips and also whether it can be used for biological control.
In the case of the Frankliniella sp. (pale morphotype), the most abundant genus was Enterobacter sp. To date, this bacterium has been reported as part of the S. dorsalis microbiota [7] but not in amounts as high as those found in this work. This bacterium has been associated with probiotic activities in the Mediterranean fly Ceratitis capitata, improving mating competitiveness, flight capacity, and longevity under starvation [73,74]. Another bacterial genus found here was Ehrlichia in the F. panamensis (dark morphotype). The presence of Ehrlichia suggests that these bacteria may be endosymbiont in pests. To date, no studies have shown the existence of this genus in thrips; however, Ehrlichia chaffeensis has been reported in other insects such as ticks, particularly in Amblyomma americanum.
Chao indices obtained from the V1-V2 region of the S. hansoni (brown morphotype) suggest that the richness of the bacterial communities is low. Furthermore, Shannon and Simpson indices were lower than Chao1, suggesting that there is community uniformity and OTU dominance in the S. hansoni samples used in this study. The evaluation of the microbiota diversity in the three morphotypes of thrips with the sequencing of the V4-V5 region of the 16S rRNA gene showed that there are no significant differences in alpha diversity indices among S. hansoni (brown morphotype), Frankliniella (pale morphotype), and F. panamensis (dark morphotype); however, regarding beta diversity and composition, the differences were very evident and strong. Similar studies have been conducted in T. tabaci suggesting that bacterial communities differ by the host [11]. Another study by Kaczmarczyk et al. [25] in H. carpathicus showed that the microbial communities varied according to the development of the life cycle, since they mostly found Wolbachia in the pupa stages.
Overall, the results obtained by culture-independent studies corroborate the findings from the culture-dependent analysis, where the presence of the Firmicutes and Proteobacteria phyla was evidenced, in addition to that of the genera Bacillus, Pantoea, Erwinia, and Sphingomonas. It is noteworthy that, due to the depth of the NGS sequencing, more genera were observed than those evidenced by culture-dependent analysis. In addition, the genus Wolbachia was found, which is an intracellular bacterium that generates symbiosis with the host and does not grow in culture medium. However, in the present study, anaerobic bacteria were not analyzed.

Conclusions
This work describes bacterial diversity present in avocado thrips species from Northeastern Antioquia, in Colombia. Our results complement the microbiota information found in various thrips species, including Frankliniella sp., F. panamensis, S. hansoni; F. occidentalis, T. tabaci, F. tritici, and H. haemorrhoidalis from avocado, bean, onion, nightshade, and wheat crops. Using molecular and traditional microbiological methods, a low bacterial diversity was found at the phylum and genus level. The application of this approach resulted in a prevalence of endosymbionts, particularly the Wolbachia, Pantoea, and Sphingomonas genera.
Additionally, this work describes the microbiota associated with three morphotypes of thrips: brown morphotype (S. hansoni), pale morphotype (Frankliniella), and dark morphotype (F. panamensis), found in avocado in Antioquia using next-generation sequencing approaches. Our results corroborate the presence of Wolbachia in the three morphotypes of thrips with prevalence in S. hansoni (brown morphotype). Furthermore, Enterobacter was found in Frankliniella sp. (pale morphotype), and Ehrlichia was found in F. panamensis (dark morphotype). Moreover, a higher α-diversity was estimated in S. hansoni.
The results obtained in this study showed that the bacterial community diversity present in adults of avocado thrips is low, having a high predominance of a few bacterial components. In addition, the microbiota in the morphotypes of thrips according to morphotype showed a large variation in the specimens analyzed. A strong predominance of the endosymbiont Wolbachia was also found here, which has been studied globally as a biological control in insects [75]. Therefore, these results are relevant to complementing the knowledge about the microbiota, as this information may be useful for biotechnological processes, including biological control methods based on symbiont employment or selective targeting.

Supplementary Materials:
The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/d14070540/s1, Figure S1: Dendrogram of Ribosomal Intergenic Spacer analysis (RISA) using complete linkage and Pearson correlation distance metric of bacterial isolates found in thrips; Figure S2: Morphology of bacterial isolates found from avocado thrips; Figure S3: Bands patterns on TGGE gels of the amplified product of the 16S rRNA gene from thrips DNA samples; Figure S4: Similarity clustering dendrogram using Dice-UPGMA of thrips samples; Figure S5: Rarefaction curve obtained for V1-V2 sequencing of the 16S rRNA gene in the brown morphotype collected from Antioquian avocado; Figure S6: Rarefaction curve obtained for V4-V5 sequencing of the 16S rRNA gene in the brown morphotype collected from Antioquian avocado; Table S1: Information on bacterial thrips samples, colony forming units (CFU), municipalities where the samples were taken and variety of avocado used for the tests; Table S2: Identification of bacteria isolated from adult thrips, according to their similarity with the sequences of the 16S RNAr gene. The identification of the species was considered positive when the % similarity was greater than 98%; Table S3: Information on used samples of thrips for independent cultivation, morphotypes, farm and avocado variety; Table S4: Taxonomic classification of bacteria and clones found in adults of thrips using temporal temperature gradient electrophoresis (TGGE) analysis, according to their similarity to the 16S rRNA gene sequences recorded in the GenBank and RDP II databases; Table S5: List of genera and percentages of relative abundance of bacteria identified in the three morphotypes with the V4-V5 region; Table S6: ANOVA test performed for the diversity indices found for the three morphotypes with V4-V5 region; Table S7: PERMAMOVA test used to compare the bacterial community between three morphotypes based on the Bray Curtis distance.

Institutional Review Board Statement:
The insect collections were made between April and August 2017 through the certificate issued by Universidad Nacional de Colombia, Sede Medellín with the permission to collect specimens of wild species of biological diversity with scientific research purposes produced by the National Authority for Environmental Licenses (ANLA) to Universidad Nacional de Colombia through resolution No. 0255 of 14 March 2014 (article 3). Thrips were collected on private property and permission was received from the landowners before sampling procedures.

Data Availability Statement:
The datasets supporting the conclusions of this article are included within the article and its additional file. The newly generated sequences are submitted to the Gen-Bank database under accession numbers MZ222250-MZ222267 for 16S rRNA sequence, MZ233678-MZ233724 TGGE band, MZ230003-MZ230010 Clones and ON012588-ON012603 for gyrase sequence. The source data of the raw Illumina sequence reads of 16S amplicons from thrips samples can be found in the Sequence Read Archive (SRA) under BioProject ID: PRJNA837112 available at http://www.ncbi.nlm.nih.gov/bioproject/837112 accessed on 30 June 2022.