Sample collection was conducted in Hoxud County in the Xinjiang Uyghur Autonomous Region of China (91°45′42″E, 40°39′75″N). The area has a semi-arid climate with a mean annual temperature of 12.56°C and mean annual precipitation of 591 mm. Two sites located in a watershed area of seasonal floods were selected for sample collection. One site is a radiation-contaminated region with ¹³⁷Cs accumulation, attributed to natural precipitation and concentration within alluvium of seasonal floods. The other site is 8 km away and not affected by radionuclide (control site). Moreover, the soil had a salt content above 2% on the surface layer (<20 cm deep). Halophyte K. schrenkianum was the most dominant plant population in the local habitat, especially at the radionuclide contaminated site.

Kalidium schrenkianum plants were collected from two sites with different radioactivity levels in April (spring), July (summer), September (autumn), and December (winter) of 2018. Plants materials were collected using a random sampling approach from a 50 m × 50 m square plot in each site. Five plants per plot were randomly selected (at least 15 m apart) and uprooted. Whole plants were labeled and placed in large autoclaved paper envelopes, then transported to the laboratory in an icebox. The samples were then stored at −80°C until further use. Additionally, five soil samples were collected in September 2018 (autumn) from the surface soil (0–20 cm deep) along the diagonal of the square plot using shovels. A total of 40 samples (2 radiation levels × 4 seasons × 2 tissue types × 5 plant replicates) were used for endophyte analysis. The soil samples from each site were sieved to remove rocks and plant litter, thoroughly mixed, labeled, and packed in cloth bags, and stored at 4°C to transport to the laboratory.

The radioactivity of ¹³⁷Cs in soil samples was analyzed at the Northwest Institute of Nuclear Technology (Xi’an, Shaanxi). One hundred grams of each soil sample were put into an environmental source box (ϕ 45 cm × 25 cm plastic box) and analyzed on an HPGe γ spectrometer (DSPEC-281, ORTEC & DSA-2000, Canberra). Each soil sample was placed 8 cm away from the detector, and the measurement time ranged from 1 to 3 days. The 661 keV peak of ¹³⁷Cs in soil samples was measured and analyzed by net counting, and the radionuclide activity of ¹³⁷Cs in the soil was calculated (Tang et al., 2008). Soil samples from the control site did not show radiation pollution (10–20 Bq/kg), while those from the radionuclide contaminated site showed significantly higher reads (>60 Bq/kg). Also, the physical and chemical indicators of the mixed soil samples from both sites were measured at the Institute of Quality Standards & Testing Technology for Agro-Products, Xinjiang Academy of Agricultural Sciences (Supplementary Table 1).

Total DNA was extracted from the aerial parts and roots of K. schrenkianum separately. Each plant material (20 g) was surface sterilized with 75% (v/v) ethanol for 1 min, 3.25% (w/v) sodium hypochlorite for 3 min, and 75% (v/v) ethanol for 30 s (Guo et al., 2000). Genomic DNA extraction was performed using the cetyltrimethylammonium bromide (CTAB) method (Guo et al., 2000). Briefly, 1 g of the surface-sterilized plant material was freeze-dried in liquid nitrogen, then ground to fine powder using a mortar and pestle. The ground samples were quickly transferred to a tube with 5 ml 2 × CTAB extraction buffer [2% (w/v), 100 mM Tris–HCl, 1.4 M NaCl, 20 mM EDTA, 1.5% (w/v) polyvinyl-pyrrolidone (PVP), 0.5% (w/v) 2-mercaptoethanol; pH 8.0; preheated to 65°C], and then incubated in a 60°C water bath for 30 min with occasional gentle swirling. Next, 500 μl of chloroform: isoamyl alcohol (24:1) was added into each tube and mixed thoroughly to form an emulsion. The mixture was centrifuged at 11,900 g for 15 min at room temperature, and the supernatant containing DNA was removed into a fresh 1.5 ml tube and re-extracted twice. Subsequently, 50 μl of 5 M KOAc was added into the supernatant followed by 400 μl of isopropanol and inverted gently to mix. The genomic DNA was precipitated overnight at 4°C and centrifuged at 9,200 g for 2 min. The DNA pellet was washed twice with 70% (v/v) ethanol and dried using SpeedVac2 (AES 1010; Savant, Holbrook, NY, United States) for 10 min or until dry. The DNA pellet was resuspended in 100 μl TE buffer (10 mM Tris–HCl, 1 mM EDTA). The DNA concentration and quality were measured using a NanoDrop 1000 Spectrophotometer (Thermo Scientific, Wilmington, NC, United States).

The V5–V7 hypervariable region of the bacterial 16S ribosomal RNA gene was amplified using 799F (AACMGGATTAGATACCCKG) (Chelius and Triplett, 2001) and 1193R (ACGTCATCCCCACCTTCC) primers (Bodenhausen et al., 2013). The fungal internal transcribed spacer region 1 (ITS1 region) of ribosomal RNA was amplified using ITS1F (CTTGGTCATTTAGAGGAAGTAA) (Gardes and Bruns, 1993) and ITS2 (GCTGCGTTCTTCATCGATGC) primers (White et al., 1990). All the PCR reactions were performed in 30 μL reactions containing 15 μL of Phusion^® High-Fidelity PCR Master Mix (New England Biolabs, Ipswich, MA, United States); 0.2 μM of forward and reverse primers, and about 10 ng template DNA. Thermal cycling conditions were as follows: initial denaturation at 98°C for 1 min, 30 cycles of denaturation at 98°C for 10 s, annealing at 50°C for 30 s, and elongation at 72°C for 30 s, then a final extension at 72°C for 5 min. The PCR products were mixed with an equal volume of 1× loading buffer (containing SYB green) and separated via 2% agarose gel electrophoresis. The expected bands were cut from the gel and purified using GeneJETTM Gel Extraction Kit (Thermo Scientific, Waltham, MA, United States). Sequencing libraries were generated using Ion Plus Fragment Library Kit 48 rxns (Thermo Scientific, Waltham, MA, United States), according to the manufacturer’s instructions. The library quality was assessed on the Qubit@ 2.0 Fluorometer (Thermo Scientific, Waltham, MA, United States). Finally, the library was sequenced on an Ion S5TM XL platform.

NGS sequencing yielded 79 samples for bacterial dataset and 78 for fungal dataset, as one sample failed to acquire sufficient reads of bacteria and two samples failed for fungi. Raw reads from the bacterial and fungal dataset were demultiplexed and quality filtered using QIIME software (Version1.7.0). Reads with a quality score <20 and those lacking complete barcode and primers were excluded from further analysis. Chimeric sequences were removed using USEARCH software. Subsequently, DADA2 workflow (Callahan et al., 2016) was used to remove singletons and doubletons. Both bacterial and fungal datasets were dereplicated to generate the amplicon sequence variants (ASVs).

The bacterial and fungal taxa were determined using the pipeline described by Zhu et al. (2021). The bacterial and fungal ASVs were first identified using Naïve Bayes approach with a minimum of 75 bootstrap calls following DADA2 workflow (Callahan et al., 2016) against SILVA version 132 (Quast et al., 2013) and UNITE general FASTA release for Fungi version 8.0 (Nilsson et al., 2019). The ASVs which were not assigned to genus (for bacteria) or species (for fungi) level were clustered into different operational taxonomic units (OTUs) based on 97% similarity. One random sequence was selected from each OTU in bacterial or fungal to assign taxonomy with SILVA or UNITE references, respectively. After assigning the taxonomy, all ASVs or OTUs that were assigned to non-bacterial, Cyanobacteria phylum, or Rickettsiales order in the bacterial dataset and to non-fungal in the fungal dataset were removed. The bacterial ASVs or OTUs were then agglomerated at the genus level, while the fungal ASVs or OTUs were agglomerated at the species level, with identical assignments using “phyloseq” package (Callahan et al., 2016).

The ASV and OTU sequences have been deposited in the GenBank of National Center for Biotechnology Information under the accession numbers MZ911967-MZ912492 (bacterial 16S sequences) and MZ919361-MZ919966 (fungal ITS1 sequences). Also, the raw sequences have been deposited in the Sequence Read Archive of NCBI under BioProject PRJNA757585 (for bacterial data) and PRJNA757584 (for fungal data).

Reviewer link:

Taxon agglomeration yielded a bacterial and fungal dataset of 2,234,413 and 3,272,864 reads, respectively. Subsequently, singletons and doubletons were filtered from the datasets, and taxa with <0.001 relative abundance were removed to reduce the noise. This resulted in a bacterial dataset of 565 taxa from 79 samples and a fungal dataset of 606 taxa from 78 samples. The bacterial and fungal datasets were then rarefied to 7,500 and 7,400 reads, respectively; the minimum value of read sums among all samples of the dataset. Finally, a bacterial dataset comprising 564 taxa and 79 samples and a fungal dataset containing 606 taxa and 78 samples were generated.

Statistical analyses were performed using R version 3.6.1 (R Development Core Team, 2016). Graphs were plotted with R packages “ggplot2” (Wickham, 2016), “grid” (Murrell, 2005), and “gridExtra” (Auguie, 2017). The distance matrices of endophytic community compositions were constructed by calculating dissimilarities using Bray–Curtis method (Faith et al., 1987). Non-metric multidimensional scaling (NMDS) was used to visualize the community composition dissimilarity of endophytic bacteria or fungi among the different seasons or radiation levels using metaMDS function in “vegan” package (Oksanen et al., 2016). Analysis of similarities (ANOSIM) was applied to analyze the differences in microbial composition between plant tissues or radiation levels. Permutational multivariate analysis of variance (PerMANOVA) with 999 permutations was implemented with adonis in “vegan” package to determine the environmental influence on microbiota composition. Community diversity of endophytic bacteria or fungi was estimated for incidence data with Hill numbers for q = 0, 1, 2 (species richness, the exponential of Shannon entropy, and the inverse of Simpson concentration) using the “iNEXT” package (Hsieh et al., 2016).

The neutral community model (NCM) and Raup-Crick index (RCI) were used to evaluate the stochasticity in the assembly of communities (Sloan et al., 2006; Chase et al., 2011). Bacterial and fungal datasets were conducted for NCM calculation and demonstration using R scripts (Burns et al., 2016) and (Chen et al., 2019). The parameter Nm in the model estimates dispersal between communities. The parameter R² represents the overall fit to the neutral model (Sloan et al., 2006). Calculation of 95% confidence intervals around all fitting statistics was done by bootstrapping with 1000 bootstrap replicates. In addition, taxa not fit to neutral model prediction were identified, as higher occupancy or abundance than expected represented plant selection or dispersal limitation (endemic), respectively (Shade and Stopnisek, 2019). RCI was calculated using the R scripts of Chase et al. (2011), and stochasticity was recognized from the proportions of community pairs that fell within | RCI| < 0.95 (Gao et al., 2020).

Correlation analysis was applied for bacterial genera, fungal species, and corresponding unassigned OTUs to deduce inter-taxa interactions. The co-occurrences between endophytic bacterial and fungal taxa in roots or aerial tissues were inspected within each treatment, namely two radiation levels and four seasons. Spearman’s rho statistic was used to estimate correlation with function cor.test in “stats” package (R Development Core Team, 2016). The frequencies of co-occurrence association (with | rho| > 0.8, p < 0.01) presented within four seasons were recorded to reveal the association persistency among seasons. If a co-occurrence association between certain taxa presented in two or more seasons, it was kept in a combined network for visualization. The combined co-occurrence networks were visualized with “igraph” package (Csardi and Nepusz, 2006). The structural dissimilarities between networks were calculated according to Schieber et al. (2017). The distance matrices were ordinated with NMDS and plotted using metaMDS function in “vegan” package as described above. Network characteristics were determined using functions in “igraph” package.

The bacterial phylum Gammaproteobacteria predominated the bacterial community in aerial tissues of K. schrenkianum in control and radiation stressed sites, followed by Alphaproteobacteria and unassigned Actinobacteria (Figure 1A). The aerial bacterial communities showed mild seasonal dynamics at the phylum level, that notable proportions of Alphaproteobacteria and unassigned Actinobacteria were observed in the control and radiation stressed sites, during autumn and summer, respectively. Similarly, the dominant phylum in the root community were Gammaproteobacteria, unassigned Actinobacteria, and Alphaproteobacteria. Like aerial communities, there was a seasonal burst of certain phyla in the control site in autumn, and the radiation stressed site in summer, including Acidimicrobiia, Nitriliruptoria, and Thermoleophilia. In the control site, their presence lasted through winter until next spring. At the bacterial genus level, Stenotrophomonas predominated aerial and root communities, followed by Ralstonia and Pseudomonas in aerial tissues and Prauserella and Actinophytocola in the roots (Figure 1B).

The fungal communities showed accordant composition in the aerial tissues and roots at the phylum level, both predominated by Dothideomycetes (Figure 2A). There were more non-Dothideomycetes fungi in the control site than in the radiation stressed site. Notably, unassigned Fungi predominated fungal communities in the roots at both sites, followed by Dothideomycetes and Sodariomycetes. At the fungal species level, Neocamarosporium spp. (F_OTU_388 and F_OTU_404) dominated the aerial fungal communities and showed dramatic colonization preference to sites (Figure 2B). The dominant taxa in the root communities were Fungi_F_OTU_767 and Fungi_F_OTU_704, which were more in the control site than the radiation stressed site. The next abundant taxa were Monosporascus_F_OTU_360 and Ascomycota_F_OTU_639, mostly inhabiting the radiation stressed site. Intriguingly, the two dominant taxa in the aerial community, Neocamarosporium_F_OTU_388 and Neocamarosporium_F_OTU_404 were also abundant in the roots, implying that the taxa are not tissue-specific. Notably, Fungi_F_OTU_767 and Fungi_F_OTU_704 might be tissue-specific because they were not observed in aerial tissues.

The bacterial community diversity in roots was higher than in aerial tissues at both sites and in all seasons (Figure 3). However, bacterial community diversity in aerial tissues at the control site was relatively low in spring and summer but sharply rose in autumn (Figure 3A). The highest diversity of bacterial communities was recorded in the roots during summer and autumn at the control site (Figure 3C). At the radiation stressed site, the highest bacterial diversity was observed in summer for both aerial and root tissues; however, the diversity of the root community dropped sharply in autumn.

The fungal community diversity showed different patterns from that of bacterial communities. In aerial tissues, the diversity was the same or higher than in the roots (Figure 4). The highest fungal diversity in aerial tissues was recorded in summer, while the lowest in winter (Figure 4A). Spring cradled the most diverse fungal communities in the roots at both sites, and aerial tissues at the radiation stressed site. However, fungal diversity levels were relatively low in the other three seasons, with a slight rise in autumn (Figures 4B–D).

Our results demonstrate that radiation and season alternation had no significant effect on the composition of endophytic bacterial communities. PerMANOVA results indicated that radiation stress did not shape bacterial communities with significant effects in aerial tissues (F₁,₃₇ = 0.58324, R² = 0.01552, p = 0.746) and roots (F₁,₃₇ = 0.59521, R² = 0.01542, p = 0.775). Season alternation significantly influenced bacterial community compositions in the roots (F₃,₃₅ = 2.2776, R² = 0.15952, p = 0.011), but not in the aerial tissues (F₃,₃₅ = 0.80269, R² = 0.06437, p = 0.651). ANOSIM results were similar to those obtained by PerMANOVA, reflecting the same patterns of environmental effects on bacterial community compositions. Radiation stress and season alternation did not affect aerial bacterial communities. However, seasonal changes affected bacterial communities in the roots (Figures 5A,C).

Compared to bacterial communities, fungal communities were more sensitive to changes in radiation levels than seasonal alternation. PerMANOVA results suggest that changes in radiation levels affected fungal communities in the aerial tissues and roots (aerial communities: F₁,₃₇ = 3.1267, R² = 0.07792, p = 0.003; root communities: F₁,₃₇ = 1.9933, R² = 0.05112, p = 0.005). Furthermore, seasonal alternation had no significant effect on fungal communities in aerial tissues (F₃,₃₅ = 1.1599, R² = 0.09043, p = 0.253) or root (F₃,₃₅ = 1.0634, R² = 0.08353, p = 0.313). ANOSIM results showed that only the site had significant effect on the fungal community compositions in the aerial tissues (Figures 5B,D). However, the fungal community compositions in the roots were significantly affected by both site and seasonal changes.

Raup-Crick index indicated generally high stochasticity in bacterial communities. The stochasticity indicator | RCI| < 0.95 accounted for more than 50% of all pairwise dissimilarities between bacterial communities in different tissue types, sites, and seasons (Figure 6A). Fewer fungal communities exhibited | RCI| < 0.95 larger than 50% (Figure 6B). RCI results suggested seasonal variation of stochasticity in endophyte communities, especially in the roots. The root bacterial communities at both sites showed lower stochasticity in the growing seasons than in winter and spring. In contrast, the highest stochasticity of fungal communities in aerial tissues or roots occurred in summer at both sites. Furthermore, the lowest stochasticity of fungal communities in the roots and aerial tissues occurred in winter and autumn, respectively.

The endophytic communities were tested by a NCM to infer their roles in the community assembly of each taxon. Endophytes with higher or lower frequencies than the prediction of neutral model were identified for each season, and the times not fit to prediction in the four seasons were recorded. Afipia and Xanthomonadaceae_B_OTU_211 showed higher frequencies in the aerial tissues than predicted in all the four seasons at control site, while Delftia and Stenotrophomonas showed higher frequencies in roots at radiation stressed site (Supplementary Table 2). Some fungal taxa, including Neocamarosporium spp. (F_OTU_388, F_OTU_404, and F_OTU_445) and Fungi spp. (F_OTU_704 and F_OTU_767) exhibited higher frequencies than predicted in all the four seasons, except in the roots at the radiation stressed site (Supplementary Table 3).

Co-occurrence networks of each endophytic community were built for each season (| rho| > 0.8, p < 0.01). Divergent distribution in structural similarities between networks was found, especially for fungal communities (Supplementary Figure 1). Combined co-occurrence networks were performed for endophytic communities of aerial tissues and roots at each site, with the co-occurrence association present in two or more seasons (Figure 7). Generally, the endophytes formed sparse networks due to low level of edge densities (Table 1). Bacterial communities in the roots and fungal communities in the aerial tissues at the control site possessed the most connected nodes and edges. The microbial communities in the roots at the control site had the highest centralized betweenness. Nevertheless, the connected endophyte taxa illustrated in network graph were all rare taxa with a relative abundance less than 0.005 (Supplementary Tables 6, 7).