To determine the relationships among symbionts, bacterial communities in S. furcifera, and environmental factors, a total of 142 S. furcifera individuals from 18 populations were collected from 14 locations across six countries (China, Thailand, Myanmar, Cambodia, Vietnam, and Laos) analyzed in this study. Each S. furcifera population comprised at least six independent samples. All S. furcifera samples were preserved in 100% alcohol and subsequently dispatched for 2bRAD-M sequencing. Information about the samples is detailed in Table 1 and Supplementary Table S1, wherein each replicate represents an independent S. furcifera adult.

Bacterial communities in S. furcifera were analyzed using the 2bRAD-M sequencing methods. 2bRAD-M sequencing and library construction were performed by Qingdao OE Biotech Co., Ltd. (Qingdao, China), following previously established protocols (Wang et al., 2012) with marginal modifications. Briefly, the DNA (100 ng) from each sample was digested using 4 U of the BcgI enzyme (NEB) at 37°C for 3 h. Next, the resulting DNA fragments were ligated with specific adaptors. A mix of 5 μL of digested DNA and 10 μL of ligation master mix, containing 0.2 μM of each adaptor and 800 U of T4 DNA ligase (NEB), underwent a ligation reaction at 4°C for 12 h.

Following ligation, the products were PCR amplified, and the resulting DNA was subjected to electrophoresis on an 8% polyacrylamide gel. An approximately 100-bp band was excised, and the DNA fragments were then diffused in DEPC water at 4°C for 12 h. A PCR test was conducted using primers bearing platform-specific barcodes to introduce sample-specific barcodes. Each PCR sample, totaling 20 μL, contained 25-ng of gel-extracted PCR product, 0.2 μM each of forward and reverse primers, 0.3 mM dNTP mix, 0.4 U Phusion high-fidelity DNA polymerase (NEB), and 1 × Phusion HF buffer. After PCR amplification and electrophoresis, the PCR products were purified using the QIAquick PCR purification kit (50) (Qiagen) and sequenced on the Illumina Nova PE150 platform.

To conduct the 2bRAD-M analysis, microbial genome data from the NCBI database, consisting of 173,165 species of fungi, bacteria, and archaea, were employed. Restriction enzymes of 16 type 2B were utilized to fragment the samples, and Perl scripts were used in this process. Thereafter, the RefSeq (GCF) number was used to assign 2bRAD-M tags with microbial genome information, including taxonomic data. Unique 2bRAD tags occurring only once in every GCF were selected as species-specific 2bRAD-M markers, forming a reference database. A detection threshold of 0.0001 (0.01%) relative abundance was set as the default (Franzosa et al., 2018).

To calculate the relative abundance of each bacterium, 2bRAD tags from all samples, post-quality control, were mapped against the 2bRAD marker database containing tags unique to the 26,163 microbial species using a built-in Perl script. To mitigate false positives in species identification, a G score was calculated for each identified species within a given sample. This score, derived from a harmonic means of read coverage of 2bRAD markers belonging to a species and the total number of possible 2bRAD markers for this species, was employed identify species while minimizing errors. A threshold G score of 5 was set to prevent false-positive microbial species discovery (Sun et al., 2022).

S: the number of reads assigned to all 2bRAD markers belonging to species i within a sample.

T: number of all 2bRAD markers of species i that have been sequenced within a sample.

The average read coverage of all 2bRAD markers for each species was computed to represent the number of individuals of a species present within a sample at a given sequencing depth. The relative abundance of a species was then calculated as the ratio of the number of microbial individuals belonging to that species against the total number of individuals from known species detected within a sample.

S: the number of reads assigned to all 2bRAD markers of species i within a sample.

T: the number of all theoretical 2bRAD markers of species i.

Additionally, five environmental factors, namely annual mean temperature, annual precipitation, longitude, latitude, and altitude, were analyzed in the study. Data for these factors were downloaded from the WorldClim website.¹ The impacts of climate and geographical factors on diversity and abundance indexes and symbiont infections were described using two structural equation models (SEM) with Satorra-Bentler correction. To limit heteroscedasticity, the log value of the “precipitation” values was used in the SEMs. SEM models were deemed acceptable when p value >0.05 and CFI > 0.95, after systematically excluding redundant pathways based on a lower AIC value. To standardize each parameter and eliminate variance, SEM coefficients were estimated through standardized transformation. Prior to analysis, data on Cardinium abundance (shown in Figure 1 left panel), Wolbachia abundance (shown in Figure 1 right panel), and microbial diversities of samples (shown in Figure 2) underwent normality testing using the Kolmogorov–Smirnov test and Levene’s test to assess the homogeneity of group variances. Data exhibiting a normal distribution were analyzed using one-way ANOVA with post-hoc Tukey HSD. In cases where the Cardinium abundance, Wolbachia abundance, and diversity data did not conform to normal distribution, they were analyzed using the Kruskal-Wallis test and Dunn’s test with Bonferroni correction for multiple comparisons. All statistical analyzes were conducted using SPSS 21.0. PCoA plots illustrating Bray-Curtis intersample distances and classification probabilities were generated using QIIME software (Caporaso et al., 2010). Pearson correlation analysis, performed using SPSS 21.0, was used to explore the relationships between different symbionts and diversity indexes. Furthermore, all graphical representations were generated using GraphPad Prism 9.0.0.

In the current study, a total of 18 S. furcifera populations were analyzed, comprising 9 populations from China and 9 populations from the South Asian region. The predominant bacterial symbionts identified in all S. furcifera populations were Wolbachia and Cardinium. Interestingly, the former was the most abundant bacteria in nearly all S. furcifera populations, being present in every S. furcifera sample. Notably, the relative abundance of Wolbachia exceeded 80% in two Chinese populations, specifically CS (85.95 ± 2.92%, Mean ± SEM) and FN (89.26 ± 3.01%, Mean ± SEM), as well as in all Laos populations, including VE (86.02 ± 1.57%, Mean ± SEM) and SA (86.02 ± 3.58%, Mean ± SEM). Nevertheless, significant variations in Wolbachia abundance were observed among different populations. Meanwhile, 3 S. furcifera populations from China (JY, YT, and TC) exhibited the lowest Wolbachia abundance, significantly lower than the Wolbachia-abundant populations such as FN, VE, and SA (Figure 1 right panel). Cardinium emerged as the second most abundant bacterium in S. furcifera (Figure 1 left panel). In several South Asia populations, including KO (63.88 ± 14.38%, Mean ± SEM), RO (59.96 ± 14.04%, Mean ± SEM), and LU (52.50 ± 15.74%, Mean ± SEM) as well as one Chinese population, YT (27.18 ± 10.38%, Mean ± SEM) Cardinium was the dominant symbiont, with its abundance in the KO population being significantly higher than that in the JY population (11.72 ± 7.78%, Mean ± SEM) (p < 0.05, Kruskal-Wallis test). Besides, Wolbachia abundance was less than 1% in the JY population, while Cardinium abundance was the highest (56.59 ± 11.90%, Mean ± SEM). Additionally, the primary symbiont Portiera, typically associated with whiteflies, was also detected in several S. furcifera populations, including TC, KO, VE, etc. (Figure 3). Principal Component Analysis was conducted to explore differences among bacterial communities of all S. furcifera populations, revealing distinct microbial communities in S. furcifera individuals from FN and JY populations compared to other samples (Supplementary Figure S1).

To elucidate variations in microbial communities among different S. furcifera populations, alpha diversity indices, including Shannon, Simpson, and Chao1, were calculated for each S. furcifera population. Considering that the three alpha indices of all populations did not follow a normal distribution, the Kruskal–Wallis test, and Dunn’s test with Bonferroni correction for multiple comparisons were adopted (Figure 2). Intriguingly, populations with low Wolbachia abundance exhibited higher alpha diversities in bacterial communities compared to those with high Wolbachia abundance. Specifically, populations with low Wolbachia abundance, such as JY, TC, and YT, displayed significantly higher Shannon and Simpson indices compared to populations with high Wolbachia abundance, including VE, SA, CX, and FN (Figure 2). The UPGMA hierarchical cluster diagram of different Sogatella furcifera populations with 2bRAD-M sequencing data was shown (Supplementary Figure S2), and the phylogenetic tree of different Sogatella furcifera populations with 2bRAD sequencing results was constructed (Supplementary Figure S3).

Pearson analysis was used to assess the correlations between symbionts and other bacteria in S. furcifera. As illustrated in Figure 4A, the presence of Wolbachia negatively impacted the abundance of various bacteria, encompassing Enterobacter, Acinetobacter, and Lysinibacillus. Notably, Wolbachia significantly and negatively influenced the abundance of Portiera in S. furcifera (r = −0.518, p < 0.05, Pearson analysis). While the abundance of Wolbachia was negatively correlated with that of Cardinium, the correlation was not statistically significant (r = −0.451, p = 0.06, Pearson analysis). Consistently, Cardinium significantly positively influenced the abundance of Pseudomonas (r = −0.804, p < 0.001, Pearson analysis) and that of various other bacteria, albeit the correlation was not statistically significant (Figure 4B).

Noteworthily, Pearson analysis exposed that the presence of Wolbachia significantly and negatively affected the three diversity indexes of S. furcifera, namely the Shannon index (r = −0.834, p < 0.001, Pearson analysis), Simpson index (r = −775, p < 0.001, Pearson analysis), and Chao1 index (r = −0.750, p < 0.001, Pearson analysis) (Figure 5). Cardinium infection also had a negative impact on the diversity indexes, although the correlations were not significant (Supplementary Figure S4).

The effects of five environmental factors (annual mean temperature, annual precipitation, latitude, longitude, and altitude) on S. furcifera microbiota were determined through Pearson analysis. While both Wolbachia and Cardinium abundance in S. furcifera increased with annual mean temperature, there was no significant correlation between temperature and the abundance of these two symbionts. Notwithstanding, temperature had a significantly negative impact on the abundance of numerous bacteria, such as Enterobacter, Lysinibacillus, and Acinetobacter (Figure 6A). Precipitation had a significant positive influence on Wolbachia abundance (r = 0.489, p < 0.05, Pearson analysis) (Figure 6B). Latitude was significantly and positively correlated with the abundance of many bacteria and was negatively correlated with that of Cardinium and Wolbachia, although the difference was not significant for the latter (Figure 6C). On the other hand, longitude was significantly and negatively correlated with the abundance of Wolbachia (r = 0.750, p < 0.001, Pearson analysis) and positively correlated with the abundance of other bacteria in S. furcifera (Figure 6D). Finally, altitude did not significantly affect bacterial abundance in S. furcifera (Supplementary Figure S5).

SEM analysis validated the significant influence of precipitation on Wolbachia abundance (r = 0.400, z = 1.967, p < 0.05, SEM) and the negative effect of latitude on Cardinium abundance (r = −3.021, z = −2.498, p < 0.05, SEM). In contrast to Pearson analysis, SEM analysis revealed a significant negative correlation between Wolbachia and Cardinium (r = −0.87, z = −2.785, p < 0.01, SEM). Lastly, temperature had a numerically positive effect on Wolbachia abundance and a negative effect on Cardinium abundance (Figure 7).

Furthermore, the effects of environmental factors on microbiota diversity in S. furcifera were explored. According to Pearson analysis, both latitude and longitude exerted a positive influence on all three diversity indexes (Shannon, Simpson, and Chao1 indexes) of bacterial communities in S. furcifera, whereas altitude had a marginal negative effect on microbiota diversity (Figure 8). In contrast, SEM analysis identified that altitude negatively impacted the Shannon index (r = −0.806, z = −2.096, p < 0.05, SEM), which differed from the results of Pearson analysis (Figure 8G). SEM analysis uncovered that temperature had a negative influence on both the Shannon index (r = −1.879, z = −2.527, p < 0.05, SEM) and Simpson index (r = −3.189, z = −3.493, p < 0.001, SEM) (Figure 9).