Bacterial communities of the psyllid pest Bactericera cockerelli (Hemiptera: Triozidae) Central haplotype of tomato crops cultivated at different locations of Mexico

Bactericera cockerelli collecting

We collected adult psyllids of B. cockerelli, sexing was not performed for this study. Psyllids were collected from potato (Solanum lycopersicum) crops under greenhouse conditions and open field cultivation in four Mexican localities between 2015–2017 (Table 1). The locations of the crops were Nuevo León (24°49′04.9″ N, 100°04′32.0″ W), Sinaloa (25°27′54.4″ N, 108°26′11.6″ W), and Coahuila (25°29′00.1″ N, 100°56′52.2″ W) in northern Mexico and Querétaro (20°34′21.7″ N, 100°25′34.2″ W) in central Mexico (Fig. 1). Psyllids in the adult stage were identified for the ground colour white yellowing in the thorax with a large well defined brown marking (Sulc, 1909; Munyaneza, 2012). The psyllid collecting was carried out during their feeding periods, placed in 95% (v/v) ethanol and stored at −20 °C until total DNA extraction.

DNA extraction

Fourteen psyllids were collected obtained two pools per localization (a total of eight samples * 7 insects × pool = 14 insects/location) (Table 1). We washed each insect pool three times with 70% (v/v) ethanol, decanting the supernatant to remove environmental contaminants before DNA extraction. After, the samples were processed as described by Caamal et al. (2019). The quality of the extracted DNA was assessed using a NanoDrop 2000c spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA). Finally, the extracted DNA was stored at −20 °C for further analysis.

Detection and identification of B. cockerelli haplotypes (mtCOI gene)

The mtCOI gene from psyllid-DNA was amplified by PCR, using primers pairs CO1F3/CO1R3 reported by Crosslin, Lin & Munyaneza (2011). PCR reactions were performed using the manufacturer’s instructions for Go Taq DNA Polymerase (Promega, Madison, WI, USA) and was used as temperate 100 ng of DNA. The thermocycling profile included an initial denaturation step for 5 min at 95 °C followed by amplification conditions for 35 cycles (95 °C, 30 s; 55 °C, 1 min, and 72 °C, 1 min), and a extension condition at 72 °C for 5 min. PCR product was purified using Zymoclean^TM Gel DNA Recovery Kit (Irvine, CA, USA), and Sanger sequenced in triplicate. The partial COI sequences obtained of Nuevo León, Sinaloa, Coahuila, and Querétaro locations were employed to identify the haplotype based on the analysis of a single nucleotide polymorphism (SNP) that existed in the mtCOI sequences between psyllids, as was described in Swisher, Munyaneza & Crosslin (2012). Previously identified sequences were used as reference: Western haplotype (JQ708095, AY971885) and Central haplotype (JQ708094, FJ175374, EF372597). Phylogenetic analyses of the sequences mtCOI gene obtained from this work were performed. The mtCOI genes sequences of psyllid used in this study were: Northwestern: JQ708093.1, KX510008.1, and KR534768.1; Central and Western: KX130767.1, JQ708094.1, JQ708095.1, FJ175374.1, EF372597.1, MK054304.1, Southwestern: KC305359.1. We used the gene sequence KP223855.1 of Bactericera maculipennis as an outgroup sequence. The phylogenetic analysis was performed with MUSCLE alignments with default parameters followed by phylogenetic tree were inferred by the maximum likelihood methodology coupled with Kimura-2 parameter model, with a bootstrapping of 1,000 replicates for statistical support of lengths of evolutionary branches, employing the Molecular Evolutionary Genetics Analysis software package version 7 (MEGA7) (Kumar, Stecher & Tamura, 2016).

Detection of Wolbachia and phylogenetic analysis of endosymbiont Wolbachia

The detection of the endosymbiont Wolbachia was carried out by PCR using the primer sequences for the gene encoding a surface protein of Wolbachia (wsp 81F/wsp 691R) reported in Zhou, Rousset & O’Neil (1998). PCR amplification was performed with the same reagents and concentrations described in the previous section, with 100 ng of DNA. The PCR program was used under the thermocycling conditions described in the previous section. The PCR products were purified using Zymoclean^TM Gel DNA Recovery Kit (Irvine, CA, USA), and Sanger sequenced in triplicate. Phylogenetic analysis of wsp was performed as described previously. The wsp gene sequences of Wolbachia used here were: Strain A: EU916190.1, GU014541.1, AB094372.1, JQ837254.1, KM267306.1, and AY971950.1; Strain B AY971930.1, AY971915.1, KM267307.1, AY971925.1, AY971917.1, AY971927.1, MK550508.1, MK534180.1, KP208732.1, GQ385975.1, KP208725.1, AB094362.1, AB094355.1, KP208729.1; and Strain c AJ252176.1, AJ252178.1, AJ252062.1; and Strain d: AJ252177.1, AJ252061.1, AJ252175.1 sequences of Wolbachia of Nematoda phylum.

16S V3 rDNA sequencing

Seven psyllids was pooled (DNA pools). A total of eight DNA pools (2 pools per location) were used to amplify the hypervariable region of the bacterial 16S rRNA. The V3 region was amplified by PCR (V3-338F and V3-533R) (Huse et al., 2008) the primers contain Illumina adapters. Sequences were obtained on the Illumina MiniSeq platform in a 2 × 150 bp (300 cycles).

Microbial diversity analysis

We applied the modified algorithm of Mott for the sequencing reads, and the threshold phred (Q) value was established at 33 (Ewing et al., 1998; Ewing & Green, 1998). Then the resulting sequences were paired-end, and we used BBmerge to merge them (Kearse et al., 2012; Bushnell, Rood & Singer, 2017). We used the Geneious assembler for de novo assembling of the merged sequences with an identity cut-off at 98% (Geneious Prime 2019.2.3; www.geneious.com). We applied Megablast for contigs, unassembled reads, and with the GenBank database non-redundant (“nt-nr”) for each pooled sample (Randle-Boggis et al., 2016). These results served as the foundation to establish a specific database on a pooled sample basis. The Sequence Classifier applied a 99% identity threshold for the taxonomical assignation of the sequence reads per pooled sample. It uses its specific database to yield the frequencies into the operational taxonomic units (OTUs) table (Geneious Prime 2019.2.3; www.geneious.com). We used the R programming language (R Core Team, 2022), starting with the use of the best hits frequencies for each OTU divided by the hits total amount for each pooled sample for data normalization and applying several packages (vegan, iNEXT, pheatmap) and tailored scripts for ecological data analyses, such as alpha indices (Chao1, Shannon, and Simpson), Welch t-test was performed to compare between the ‘Ca. Liberibacter’ relative abundance (high- and low-carrier), Bray-Curtis dissimilarity distance estimations, and principal coordinates analysis were performed for beta diversity analyzes. The PERMANOVA statistical analysis was performed to determine the significance of microbial community differences. And, the differential abundance OTUs (DAOTUs) analysis of the pooled samples (Oksanen et al., 2014; Hsieh, Ma & Chao, 2016).

Sequence analysis of mtCOI gene for identification of haplotype of B. cockerelli

All psyllids’ samples were assessed through PCR using the primers specific to the B. cockerelli mtCOI gene. The resulting PCR products for all samples were positive with an amplicon length of approximately 455 bp (Fig. 2A). The mtCOI amplicon was sequenced and subjected to analysis to identify the haplotype of the psyllids in search of only one single-nucleotide polymorphism (SNP) at nucleotide 51 that exists between the Western (nucleotide 51:C) and Central (nucleotide 51:T) haplotypes in a 455-bp fragment of the mtCOI gene (Fig. S1). With the consensus sequence length of 453 nucleotides, the mtCOI sequences obtained in this work align from the nucleotide 48 up to nucleotide 500 with the reference mtCOI sequences numbers: Western (JQ708095 and AY971885), and Central (JQ708094, FJ175374, and EF372597). All the samples of this study had thymine in the residue corresponding to the nucleotide 51 of the reference sequences represented the Central haplotype of the psyllid (Fig. S1).

In the maximum likelihood (ML) tree, mtCOI sequences identified in the present study formed a robustly supported clade (bootstrap: 94%) with those of psyllids belonging to the Central and Western haplotype, with a close relationship with the Central haplotype (JQ708094.1 and MK054304.1), consistent with the results of the analysis of SNP (Fig. 2B).

Detection and identification of symbionts Wolbachia from B. cockerelli Central haplotype

The Alphaproteobacteria Wolbachia was detected using PCR that amplifies a fragment of the gene encoding a surface protein of Wolbachia (wsp) in all B. cockerelli samples collected in this work (Fig. 3A). The PCR products were sequenced, and the wsp sequences were analyzed with blastn algorithm for global alignment, with 100% identity to Wolbachia endosymbiont of B. cockerelli (GenBank accession number: KM267307; Cooper et al., 2015) (Fig. 3A). The phylogenetic analysis clustered the Wolbachia of B. cockerelli (Central haplotype) from Querétaro, Coahuila, Nuevo León, and Sinaloa with the monophyletic grouping of Wolbachia strain B from B. cockerelli (GenBank accession number: KM267307; Cooper et al., 2015) (Fig. 3B). In addition, the sequences from the samples analyzed here were related to the monophyletic grouping Wolbachia strain B from other insects (e.g., GenBank accession number: GQ385975.1, AB094355, and AB094362.1), and the paraphyletic group of Wolbachia strain A identified in B. cockerelli (GenBank accession number: KM267306; Cooper et al., 2015) (Fig. 3B).

Sequencing run metrics and microbial community diversity analyzes from B. cockerelli

We obtained 830,240 reads of eight samples of adult psyllids collected in Querétaro, Coahuila, Nuevo León, and Sinaloa. After bioinformatic processing (trimming, quality control, and paired-end merging) the whole sequences we kept was 399,936 reads that, in turn, were classified into 689 OTUs at 99% sequence identity (Table 2). It is worth noting that B. cockerelli samples might carry bacterial plant pathogens. Among the identified OTUs, we found the presence of the bacterial plant pathogen ‘Candidatus Liberibacter’. Moreover, the OTU corresponding to ‘Candidatus Liberibacter’ is among the most abundant OTUs. We performed an unsupervised hierarchical biclustering analysis of the most abundant OTUs (>1% relative abundance/sample). This analysis showed the grouping of the samples mainly based on the relative abundance of eight OTUs (high and low relative abundance), among which ‘Candidatus Liberibacter’ OTU is present (Fig. 4). Then, we assigned the samples into two well-defined groups, high- (BcAdOFCoa01, BcAdOFCoa02, BcAdGHNL02, and BcAdGHQue01) and low-carriers (BcAdGHNL01, BcAdGHQue02, BcAdOFSin01, and BcAdOFSin02), based upon the relative abundance of ‘Candidatus Liberibacter’. Moreover, we assessed the alpha diversity indices through inter- and extrapolation sampling curves analysis. The estimation of diversity (Chao1 index, q = 0), richness (Shannon index, q = 1), and evenness (Simpson index, q = 2) showed differences among insect samples also due to ‘Candidatus Liberibacter’ relative abundance (high- and low-carrier) independently of the location (Fig. 5). Also, the overall alpha diversity indices were estimated from the bacterial community abundance to assess diversity (observed OTUs and Chao1 index), richness (Shannon index), and evenness (Simpson index) for the B. cockerelli Central haplotype samples. For all samples, the observed OTUs and Chao1 index, Shannon index, and Simpson index indicated a high specialization of the bacterial community associated with B. cockerelli (85 ± 16.80, 165.96 ± 29.935, 1.328 ± 0.472, and 0.582 ± 0.167, respectively) (Table 2 and Fig. S2). It was also possible to infer a grouping based on ‘Candidatus Liberibacter’ relative abundance since the Shannon index values above 1.417 to 1.951 and the Simpson index values above 0.6362 to 0.7996 corresponding to the samples with high relative abundance of ‘Candidatus Liberibacter’ (Table 2 and Fig. S2).

We proceeded with the sample clustering analysis or β-diversity analysis applying a Principal Coordinate analysis (PCoA) based on a Bray-Curtis dissimilarity distances estimations between samples. The PCoA analysis showed two discrete groups based on the relative abundance of ‘Candidatus Liberibacter’ (high- and low-carriers) from B. cockerelli Central haplotype samples (Fig. 6). Furthermore, we carried out the Permutational multivariate analysis of variance (PERMANOVA) to determine whether the relative abundance of ‘Candidatus Liberibacter’ (high- and low-carriers) or location and even crop conditions (open field or greenhouse) played a significant role in the clustering of the samples. The PERMANOVA results showed that ‘Candidatus Liberibacter’ relative abundance (high- and low-carriers) was the main factor influencing the bacterial communities’ assemblages of B. cockerelli central haplotypes samples (R² = 0.58536, P = 0.025), and location and crop condition (greenhouse or open field) did not have any significant effect (R² = 0.45264, P = 0.39; and R² = 0.03368, P = 1, respectively). Furthermore, we extended the PERMANOVA analysis to determine the interactions between ‘Candidatus Liberibacter’ relative abundance and location and ‘Candidatus Liberibacter’ relative abundance and crop condition. Both PERMANOVA interactions analyzes showed that ‘Candidatus Liberibacter’ relative abundance (relative abundance:location, R² = 0.58536, P = 0.048; relative abundance:crop condition, R² = 0.58536, P = 0.019) was the main factor that affected the B. cockerelli microbial structure assemblage above location and crop condition (R² = 0.15666, P = 0.844; and R² = 0.03680, P = 0.665, respectively).

The microbial community structures of B. cockerelli central haplotype across all samples showed that the dominant phylum was Proteobacteria with relative abundances ranging from 91 to 99%. At Class level the dominant class was Alphaproteobacteria with a relative abundance ranging from 95 to 99%. At the genus level the dominance of the Proteobacteria phylum followed the same trend (relative abundance from 87 to 99%), and the most abundant genera were Neorickettsia and ‘Candidatus Liberibacter’ ranging from 18% to 98%, and from 0.01% to 66%. The ‘Candidatus Liberibacter’ high-carrier samples (BcAdOFCoa01, BcAdOFCoa02, BcAdGHNL02, and BcAdGHQue01) showed relative abundances ranging from 51% to 66%, giving rise to be the most abundant genus in those samples, and the low-carrier samples (BcAdGHNL01, BcAdGHQue02, BcAdOFSin01, and BcAdOFSin02) showed relative abundances ranging from 0.01% to 2.26%, in these samples Neorickettsia genus was dominant with relative abundances ranging from 64 to 98% (Fig. 7A). To determine which OTUs underwent significant relative abundance changes in the microbial assemblage across all samples, we carried out a differential abundance OTUs (DAOTUs) analysis. Surprisingly, from the 689 OTUs assessed through the DAOTUs analysis, only the OTU corresponding to ´Candidatus Liberibacter´ showed significant differential abundance between high- and low-carrier samples (log2foldchange = 6.10, P value = 0.000127, Padj = 0.0409) (Fig. 7B). Altogether, the presence and the relative abundance of ‘Candidatus Liberibacter’ suggest a pivotal role in microbial assemblage modelling.