Citrus root samples were collected from 28 distinct citrus groves situated across nine citrus-producing countries (Supplementary Table S1), encompassing all six continents that citrus can grow. The root samples were collected from local citrus varieties (scion/rootstock) in each location (Supplementary Table S1). For each country, the most representative citrus growing locations were selected for samples collection, except for Oman, which had five locations. In each grove, four healthy trees were selected and, from each of these trees, fine roots (approximately 1 mm diameter) from a depth of 5–15 cm were collected. For each tree, the samples were collected from four ordinate directions approximately 1 meter away from the trunk. The roots were removed from the soil with a shovel and then gently shaken to remove the soil that was not tightly attached to the roots. The roots from the four locations/tree were pooled and washed thrice using PBS buffer. 2 grams of washed roots were subjected for DNA extraction using a MoBio Powersoil DNA extraction kit (MoBio Laboratories Inc. Carlsbad, CA, USA) according to the manufacturer’s instructions. The DNA quality and quantity were tested using a NanoDrop device (Thermo Scientific, Wilmington, DE, USA) and through electrophoresis (0.8% agarose gel, including a 1 kb plus ladder). The DNA samples from the four trees collected from the same grove were pooled together and stored at −80°C until use. It is noteworthy that the root samples (including rhizoplane and endosphere) were collected with the rhizosphere (soil surrounding the roots), and bulk soil (Supplementary Table S1), which allows comparison with rhizosphere and bulk soil samples from same locations. Analyses of the bulk soil and rhizosphere samples were already reported in our previous study (Xu et al., 2018).

To investigate both prokaryotic and eukaryotic microbes in the root samples, 16S and ITS2 amplicon library preparation and sequencing were performed according to the manufacturer’s protocol at BGI-Shenzhen. For targeted metabarcoding sequencing, DNA fragments were amplified using primers targeting the 16S rDNA V4 region (515F and 806R) (Caporaso et al., 2011) and fungal ITS2 (Tedersoo et al., 2015). After quality control, quantification, and normalization of the DNA libraries, 250 bp paired-end reads were generated using the Illumina MiSeq250 sequencing platform according to the manufacturer’s instructions. The raw sequencing reads of root samples were deposited in the NCBI Short Read Archive BioProject PRJNA844917. For the comparative analysis with rhizosphere and soil, raw sequencing reads of the corresponding rhizosphere and bulk soil samples were obtained from Xu et al. (2018) (BioProject PRJNA362455).

The 16S rRNA gene and ITS2 amplicon sequencing data of the root samples were analyzed using the DADA2 pipeline (R package dada2, v.1.8.0) (Callahan et al., 2016) in R software (v.4.2) as follows. To compare with corresponding rhizosphere and bulk soil samples from previous study (Xu et al., 2018), we also included the bacterial 16S rDNA and fungal ITS2 amplicon data of rhizosphere and bulk soil samples for further analyses using the same methods as root samples. In detail, the raw sequence reads were first filtered, de-replicated and de-noised using DADA2 recommended parameters. Paired-end sequences were then merged, and chimeras were removed. To obtain the taxonomic information, representative sequences of ASVs were searched against SILVA database for 16S and UNITE database for ITS2 using the RDP Naive Bayesian Classifier algorithm (Wang et al., 2007). ASVs assigned with no kingdom-level classification or defined as “Unknown” at the phylum rank were removed. ASVs of 16S data classified as “Chloroplast,” “Mitochondria” or “Eukaryota” and ASVs of ITS data classified as “Bacteria” or “Plant” were further removed. Sequence alignment was performed with AlignSeqs function from DECIPHER package (v.2.20.0). The phylogenetic tree was built using phangorn package (v.2.8.1) and then rooted for downstream reproducible analysis. The relative abundance tables for taxa were generated based on the read count for each taxon across samples using the total-sum scaling (TSS) method (Weiss et al., 2017). Diversity indices were calculated using phyloseq package (v.1.40.0) (McMurdie and Holmes, 2013) and visualized using ggplot2 package (v.3.3.5). Specifically, alpha diversity indices were calculated for each sample with not-rarefied data (McMurdie and Holmes, 2014). Beta diversity was performed using Principal Coordinates Analysis (PCoA) method with Unweighted distance in not-rarefied data, considering phylogenetic relatedness between samples in calculations.

Significant differences across locations were determined with one-way Analysis of Variance (ANOVA) followed by Tukey HSD post hoc test at p-value<0.05, prior checked the normal distribution of variable by Shapiro–Wilk test. Samples related to root and corresponding rhizosphere associated to different rootstocks were grouped according to the shared rootstocks (Supplementary Table S1). Specifically, six rootstocks which were present in at least three locations were selected: Swingle citrumelo [Citrus paradisi Macf. × P. trifoliata (L.) Raf.], Rangpur lime [C. limonia L. Osb.], Trifoliate orange [C. trifoliata (L.) Raf.], Ortanique tangor [C. sinensis L. Osb. × C. reticulata Bl.], Citrange [C. sinensis (L.) Osb. × C. trifoliata (L.) Raf.] ‘Carrizo’ and ‘Troyer’. Samples associated with other rootstocks with less than three replicates were discarded from further analyses. The significant differences in alpha diversity across compartments of same locations and rootstock genotypes were determined by one-way ANOVA and Fisher’s least significant difference (LSD) post hoc test at p-value<0.05. In addition, two-way ANOVA effect with interaction between the fixed variables compartment, location and rootstock genotypes was performed. Statistical analyses of beta diversity were performed by permutational multivariate analysis of variance (PERMANOVA), using a pairwise multilevel comparison throughout Adonis (Anderson, 2001) command in vegan package (v.2.5.7), with a permutation number of 999 available.

Based on the read count abundance profiles, the features (phylum, genus) with significantly differential abundances across compartments (root versus corresponding rhizosphere) were determined using DESeq2 tool (v.1.34.0) (Love et al., 2014), fitting the negative binomial model and testing the significance with Wald test. Analyses of Compositions of Microbiome with Bias Correction (ANCOM-BC) among rootstocks for both ITS and 16S, in Root and Rhizosphere at Class level were performed usingANCOMBC package (Lin and Peddada, 2020). p-values for multiple tests were corrected using the BH method. All items with adjusted p-value<0.05 were considered significant. Moreover, based on the relative abundance of bacterial and fungi ASVs, the pairwise spearman’s correlation coefficient (ρ) values were calculated. Only correlations between two ASVs were considered as statistically robust if |ρ| ≥0.65 (Spearman’s rank correlation coefficient) and the p-value<0.05. The degree, betweenness, and modularity of each network were calculated using igraph package (v. 1.5.1). The co-occurrence networks were visualized using the ggnet2 package (v. 0.1.0).

Core members of the root and rhizosphere were calculated using the occurrence of association across samples as the criterion (Neu et al., 2021) using Core function in microbiome package (v.1.5.28) (Lahti and Shetty, 2017). Specifically, the genera that were present in more than 75% of samples across the globe for each compartment were considered as core genera (Xu et al., 2018). To generate the potential functions of core root microbiome, we first identified the genes that belonged to these core genera from our global citrus rhizosphere microbiome gene sets (Xu et al., 2018). Then according to the functional annotation of identified genes, the potential functions and KEGG pathways of core genera of root microbiome were inferred. The enrichment analysis of functional pathways compared to the functions of whole gene set were performed using fishers’ exact test based on p-value was less than 0.05.

The clean dataset obtained from root samples comprised a total of 28 samples from nine countries spanning in all the six continents (Supplementary Table S1) where citrus grows for 16S rDNA amplicon analysis and 27 samples for ITS2 amplicon analysis. A sample with fewer than 1,000 high-quality reads was excluded from the ITS dataset. The dataset contained an average of 8,166 and 17,930 read counts per sample for 16S and ITS, respectively (Supplementary Tables S2, S3). A total of 2,586 distinct prokaryotic (bacterial and archaeal) ASVs encompassing 31 phyla and 2,812 distinct fungal ASVs encompassing 12 phyla were identified using DADA2 pipeline (Supplementary Table S4). We defined the prokaryotic and fungal phyla that exhibited relative abundance >1% in at least one location as the predominant ones (prokaryotic, Figure 1A; and fungal, Figure 1B). Ten predominant prokaryotic phyla, including 9 bacterial phyla and 1 archaeal phylum (Thaumarchaeota), were identified in the citrus root samples, with Proteobacteria, Actinobacteria, Acidobacteria, and Bacteroidetes being the top dominant phyla identified (Figure 1A). Regarding the fungal phyla, the representatives included Ascomycota, Basidiomycota and, to a lesser extent, Glomeromycota and Mortierellomycota, together accounting for more than 90% of the total phyla in each location (Figure 1B).

To gain insights in the assembly clue of root microbiota, we analyzed the 16S rDNA and ITS2 amplicon sequencing data from the corresponding rhizosphere and bulk soil samples from the same trees (Xu et al., 2018) using the identical analysis pipeline as mentioned above. The alpha-diversity of both the root prokaryotic and fungal communities was lower compared to the adjacent rhizosphere samples. This was evident from the results obtained through both location-specific and overall comparative analyses (Figures 1C,D; Supplementary Figure S1, and Supplementary Tables S5, S6) (p-value <0.05; One-way ANOVA followed by Fisher’s least significant difference test). Beta-diversity analysis revealed that the root samples clearly clustered distinctly apart from the bulk soil and rhizosphere, forming two separate groups for both prokaryotic (Figure 1E, R² = 0.05; p < 0.005, PERMANOVA) and fungal (Figure 1F, R² = 0.09, p < 0.001, PERMANOVA) communities. Moreover, network analyses were performed to compare the inter-member interaction attributes between the root and rhizosphere microbial communities (Supplementary Datas S1–S4). The root prokaryotic network consisted of 334 nodes with an average degree number of 5.84 and modularity score of 0.697, which suggested that the root prokaryotic community formed a modular structure based on the modularity score threshold of 0.634, as described by Newman (2006); the corresponding rhizosphere co-network consisted of 781 nodes with an average degree number of 14.18 and modularity score of 0.567, indicative of a non-modular structure (Figures 2A,C). The root fungal network comprised of 197 nodes, with an average degree number of 4.72 and modularity score of 0.766, while the rhizosphere network comprised of 385 nodes, with an average degree number of 10.62 and modularity score of 0.574 (Figures 2B,D). These results suggested that the root prokaryotic and fungal communities were less diverse and complex but more modular than the rhizosphere ones.

In the rhizosphere and bulk soil microbiota, more predominant prokaryotic (12 phyla for rhizosphere and 13 phyla for the bulk soil) and fungal (five phyla for rhizosphere and six phyla for the bulk soil) phyla were identified than the root microbiota (Figures 1A,B). Further comparative microbiota analyses were conducted, and the results suggested that no bacterial phylum was enriched, and 20 phyla were depleted in the root microbiota compared with the corresponding rhizosphere (corrected p-value<0.05, DESeq2, Supplementary Table S7). One fungal phylum, Glomeromycota, was enriched in the root niche, while Kickxellomycota, Mortierellomycota, Chytridiomycota, Monoblepharomycota, Aphelidiomycota and Blastocladiomycota were depleted in the root microbiota compared with the rhizosphere microbiota (corrected p-value<0.05, DESeq2, Supplementary Table S7). At the genus level, 30 root-enriched and 150 root-depleted genera in bacterial community were identified (corrected p-value<0.05, DESeq2, Supplementary Data 5). For the fungal community, 21 genera were enriched, and 147 genera were depleted in the root niche compared with the rhizosphere (corrected p-value<0.05, DESeq2, Supplementary Data 5). Moreover, 62 bacterial genera and 61 fungal genera were specifically identified in the roots but not identified in the rhizosphere (Supplementary Table S8).

Interestingly, while the alpha-diversity of the prokaryotic community of citrus rhizosphere (F_7,22 = 0.60; p = 0.746) and bulk soil (F_6,18 = 0.56, p = 0.751) exhibited minimal variations across different geographic locations, the alpha-diversity of the root prokaryotic community displayed significant fluctuations across these geographic areas (F_8,27 = 6.73, p < 0.05) (Figure 1C and Supplementary Table S5). For the fungal communities, no statistically significant differences were identified for the root and bulk soil samples across all locations (Figure 1D) (F_8,25 = 2.32; p = 0.069, root; F_5,13 = 2.53, p = 0.117, bulk soil); however, the rhizosphere fungal community exhibited significant variations in term of alpha-diversity (F_7,20 = 3.80; p < 0.05) (Figure 1D and Supplementary Table S6). Because citrus genotype is one of the most significant variables in different locations, we further conducted citrus genotype-based group-wise comparisons. The results revealed that the citrus genotype indeed had strong influences on the root prokaryotic and fungal communities. The alpha diversity of root prokaryotic (F_5,21 = 9.41, p < 0.05) and fungal (F_5,20 = 3.05, p < 0.05) communities were significantly different among the six citrus genotypes (Figure 3A), while no significant difference was observed for the rhizosphere prokaryotic (F_5,21 = 1.21, p > 0.05) and fungal (F_5,19 = 0.903, p > 0.05) communities (Supplementary Figure S2A). Moreover, the beta-diversity comparison demonstrated that the factor “genotype” exerted a pronounced influence on the root prokaryotic (R² = 0.36, p-value<0.001, PERMANOVA) and fungal community (R² = 0.35, p-value<0.001, PERMANOVA).

Several predominant taxa exhibited dramatically varied relative abundance among different citrus genotypes (Figure 3B; Supplementary Figure S2B, and Supplementary Table S9). For instance, the relative abundance of class Gammaproteobacteria and Bacteroidia, two of the dominant classes in both root and rhizosphere prokaryotic community, were higher in Carrizo citrange than Trifoliate orange, while the relative abundance of class Acidobacteriia exhibited an opposite trend with relative abundance higher in Trifoliate orange than Carrizo citrange (Figure 3B; Supplementary Figure S2B, and Supplementary Table S9). Similarly, variations in the relative abundance of several taxa in the root and rhizosphere fungal community were also observed (Figure 3B; Supplementary Figure S2B, and Supplementary Table S9).

In this study, the taxa that were present in >75% of the root samples were defined as core taxa of citrus root microbiota, and 12 core prokaryotic genera and 7 core fungal genera were identified (Supplementary Data 6). The core prokaryotic genera included Acidibacter, Allorhizobium, Bradyrhizobium, Chitinophaga, Cupriavidus, Devosia, Dongia, Niastella, Pseudomonas, Sphingobium, Steroidobacter, and Streptomyces. All of these genera were among the most abundant genera in root samples (Figure 4A). Moreover, all the identified core genera, which contained several ASVs, were present within the co-occurrence network, and particularly, the ASVs members belonging to genus Pseudomonas were identified as having both high number of degree (number of connections of each node) and high values of betweenness centrality (importance for network connectivity), representing putative hubs and key connectors (Supplementary Data 1). Eleven of the 12 genera were also identified as the core members in the citrus rhizosphere microbiota (Supplementary Data 7), and Allorhizobium exhibited significantly higher relative abundance in the root samples than in the rhizosphere ones (corrected p-value<0.05, DESeq2; Supplementary Data 1). Of note, Chitinophaga was identified as a core genus in the root microbiota and not in the rhizosphere microbiota. Furthermore, Chitinophaga only exhibited relatively low relative abundance in 17 of 23 rhizosphere samples (Supplementary Datas 6, 7). The core fungal genera in the citrus root microbiota included Acrocalymma, Cladosporium, Fusarium, Gibberella, Mortierella, Neocosmospora and Volutella (Supplementary Data 6). These fungal genera were also among the most abundant fungal genera identified in the root compartment across globe (Figure 4B), and all identified taxa were present within the root co-occurrence networks and majority of them exhibited high scores of betweenness centrality (Supplementary Data 2). All fungal genera were also identified as core genera in the citrus rhizosphere microbiota. Among them, Neocosmospora exhibited significantly higher relative abundance in the root samples than in the rhizosphere ones (corrected p-value<0.05, DESeq2; Supplementary Data 5); and Fusarium, Gibberella and Mortierella exhibited relatively lower abundance in the roots (Supplementary Datas 6, 7).

Next, the potential functions of these core genera of root microbiota were predicted based on our previously established gene catalog of the global citrus rhizosphere microbiome (Xu et al., 2018). The predicted functions were distributed in 211 level 3 KEGG pathways (Supplementary Data 8). The functions of these root core genera were enriched in several KEGG pathways, including “ABC transporters,” “quorum sensing,” “biofilm formation,” “bacterial motility proteins,” “arginine and proline metabolism,” and “tyrosine metabolism,” compared with the rhizosphere (p < 0.05, fishers’ exact test) (Figure 5). On the other hand, the pathways “peptidases,” “amino sugar and nucleotide sugar metabolism,” and “starch and sucrose metabolism” were depleted in the root samples compared with the rhizosphere ones (p < 0.05, fishers’ exact test) (Figure 5).