Two mulberry cultivars, “Da Shi” (DS) and “Xuan792” (X792), were used in the study. X792 is resistant to low temperature (frozen shoot rate 5.06%), while DS is sensitive to low temperature (frozen shoot rate 30.43%). Both varieties were sampled at the experimental farm of Shandong Sericulture Research Institute (37°08′25.44′N, 121°08′33.98′E) in Yantai, Shandong Province, China, and shared the same climatic conditions. The local area has a temperate maritime climate, possessing a mean air temperature of 14°C and average annual rainfall of approximately 700 mm, with a frost-free period of more than 210 days, and an annual sunshine duration of more than 2,100 h. The maximum and the minimum temperatures in January were 5.55 and −1.97°C, respectively, and in February were 6.97 and −0.10°C, respectively. The temperature of the sampling day in January was −4~3°C, and in February was 3~9°C.

Three replicates for each sample (2 cultivars × 2 months× 2 organs) were performed in the study. The samples were collected from 10-year-old trees of mulberry cultivars, with three healthy-looking, medium-growth plants randomly selected from each site in January and February of 2019. Plants were at least 2 m apart within an area of 100 m². Healthy branches approximately 180.0 cm in length and 1.5–1.8 cm in diameter, and roots approximately 30.0 cm in length and 0.3–0.8 cm in diameter, were collected from the three individual plants of each cultivar. Whole branches and roots were placed in Ziploc bags and stored at 4°C during transportation to the laboratory, then were processed within 24 h of collection. Branches and roots were washed in running tap water to remove surface debris and then the middle 60 cm of the branches were cut into several 5.0-cm segments, while the middle 18 cm of the roots were cut into several 2.0-cm segments (15–20 cm depth). Six surface-sterilized segments of each cultivar were randomly selected, pooled, and served as one replicate for further endophyte enrichment (Nxumalo et al., 2020).

Plant materials were surface sterilized using the procedure of Du et al. (2020). Briefly, the stem and root were washed thoroughly with sterile water, then immersed in 70% ethanol for 3 min, washed with fresh sodium hypochlorite solution (2.5% available Cl⁻) for 5 min with agitation, rinsed three times with 70% ethanol for 30 s, and finally washed five times with sterile distilled water. The sterile distilled water used in the final wash was cultivated to determine the success of the surface disinfection. Briefly, 100 μl of the final rinse water was plated on an LB medium and examined for bacterial growth after incubation at 30°C for 72 h. If there was no bacterial growth, the surface-sterilization procedure was confirmed to be effective and the samples were used for further analysis.

Genomic DNA of the microbial community of the mulberry plants was extracted from stems and roots using the E.Z.N.A.® soil DNA Kit (Omega Bio-tek, Norcross, GA, United States) according to the manufacturer’s instructions. The quality of the extracted DNA was checked on a 1% agarose gel, and DNA concentration and purity were determined with a NanoDrop 2000 UV–vis spectrophotometer (Thermo Scientific, Wilmington, DE, United States). The hypervariable region V5–V7 of bacterial 16S rRNA genes was amplified by PCR using primer pairs 799F (5′-AACMGGATTAGATACCCKG-3′) and 1193R (5′-ACGTCATCCCCACCTTCC-3′; Beckers et al., 2017). The PCR mixtures contained 4 μl 5× TransStartFastPfu buffer, 2 μl of 2.5 mM dNTPs, 0.8 μl each primer (5 μM each), 0.4 μl TransStartFastPfu DNA Polymerase, 10 ng template DNA, and ddH₂O to 20 μl, and reactions were performed in triplicate. PCR cycling conditions comprised an initial denaturation at 95°C for 3 min, 27 cycles of denaturing at 95°C for 30 s, annealing at 55°C for 30 s, and extension at 72°C for 45 s, followed by a single extension at 72°C for 10 min and a continued hold at 4°C. The complete sequences generated in this study are available in the NCBI SRA database under accession number SRR18790631–SRR18790654.

The resulting PCR products were extracted from 2% agarose gels, purified using the AxyPrep DNA Gel Extraction Kit (Axygen Biosciences, Union City, CA, United States) according to the manufacturer’s instructions, and quantified using a Quantus™ Fluorometer (Promega, United States). Purified amplicons were pooled in equimolar ratios and paired-end sequenced by Majorbio Bio-Pharm Technology Co. Ltd. (Shanghai, China) using an Illumina MiSeq PE300 platform (Illumina, San Diego, CA, United States) according to standard protocols.

Raw 16S rRNA gene sequencing reads were demultiplexed, quality-filtered by fastq version 0.20.0 (Chen et al., 2018), and merged by FLASH version 1.2.7 (Magoč and Salzberg, 2011). Operational taxonomic units (OTUs) with a 97% similarity cut-off (Stackebrandt and Goebel, 1994; Edgar, 2013) were clustered using UPARSE version 7.1 (Edgar, 2013), and chimeric sequences were identified and removed. The taxonomy of each OTU representative sequence was analyzed by RDP Classifier version 2.2 (Wang et al., 2007) against the 16S rRNA database using a confidence threshold of 0.7 (Beckers et al., 2016).

Bacterial relative abundance, α-diversity, community composition, β-diversity, network structure, and functional analysis were performed using the free online Majorbio Cloud Platform.¹ Alpha diversity was calculated including Chao, ACE, Shannon, and Simpson indices using Mothur software (versionv.1.30.2). Rarefaction curves were also generated using Mothur at a 97% identity level. The Venn, Bar, and Heatmap diagram was generated using R script (version 3.3.1), and Silva (Release138; http://www.arb-silva.de) was used for taxonomic classification. β-diversity was visualized using principal coordinates analysis (PCoA) and non-metric multidimensional scaling (NMDS) based on the distance matrix, with the calculation of the Euclidean and Bray-Curtis algorithm, respectively. The potential function of bacterial communities was predicted using PICRUSt2 (Douglas et al., 2020; Phylogenetic Investigation of Communities by Reconstruction of Unobserved States), and the histogram was created using Graphypad Prism 9.3.1. The OTU abundance table was first normalised by PICRUSt2 (the PICRUSt process stores the COG information and KO information corresponding to the OTU), i.e., the influence of the number of copies of the 16S marker gene in the species genome was removed; then obtain the COG family information and KEGG Ortholog (KO) information corresponding to the OTU, and the abundance of each COG and the KO abundance were calculated; According to the information in the COG database, the description information of each COG and its function information can be parsed from the eggNOG database, so as to obtain the functional abundance spectrum; based on the information in the KEGG database, the KO, Pathway, and EC information can be obtained, and the abundance of each functional class can be calculated based on the OTU abundance. In addition, for Pathway, three levels of metabolic pathway information can be obtained by using PICRUSt2, and the abundance table of each level can be obtained respectively, significant differences were analyzed by TTEST in excel. The community differences for the 10 most abundant bacterial genera distributions were evaluated using the Student’s T and one-way ANOVA, with values of p < 0.05 considered statistically significant. Finally, a network analysis was performed by spearman using Network × software to explore the complexity of the interactions among the microbial taxa, the absolute value of the correlation coefficient ≥0.7, p < 0.05, and the diagram was drawn by Cytoscape 3.5.1. All experiments were conducted in triplicate. Data are presented as means with SDs. Statistical analyses were performed using DPS Statistics 18.10 software.² Differences between the means of different treatments were determined using the Duncan test at p < 0.05. The OTU table with taxonomic annotations was provided as Supplementary Table S1.

After read-quality filtering, a total of 1,059,743 high-quality sequences remained and were queried. The total number of bases was 399,174,839, and the average read length was 376.63 bp, (Supplementary Table S2). Rarefaction curves (Supplementary Figure S1), combined with the estimated coverage values (Table 1), suggested that the data were sufficiently large to capture most of the bacterial diversity in the samples. The number of OTUs obtained was highest in the stem of X792 in February, followed by that in the root of X792 in January, while the lowest number of OTUs was present in the stem of X792 in January.

The number of common and unique bacterial OTUs in the different samples was presented in Venn diagrams (Figure 1). The numbers of shared OTUs between January and February samples (983, Figure 1A) and between DS and X792 samples (932, Figure 1B) were both higher than the number of shared OTUs between stem and root samples (801, Figure 1C). For unique OTUs, the number obtained in January (332) was lower than the number obtained in February (667; Figure 1A); the number in the sensitive mulberry cultivar DS (302) was lower than the number in the resistant mulberry cultivar X792 (748; Figure 1B); and the number obtained from the stems (821) was higher than the number obtained from the roots (360; Figure 1C).

In addition, the number of OTUs shared between different organs and months for the DS cultivar (103; Figure 1D) was lower than that of the X792 cultivar (176; Figure 1E); the number shared between different organs and cultivars in January (97; Figure 1F) was lower compared with that in February (163; Figure 1G); and the number of OTUs shared between different months and cultivars for the stem samples (172; Figure 1H) was lower than that obtained for the root samples (260; Figure 1I). The number of unique OTUs in the roots of the DS cultivar in February (FRDS, 88) was lower than those of the other samples of the DS cultivar (JSDS, 178; JRDS, 145; and FSDS,174; Figure 1D); but when compared within stems (Figure 1H), or within February (Figure 1G), or within X792 (Figure 1E), the number of unique OTUs in the stems of X792 in February was the highest (516 in Figure 1E; 548 in Figure 1G; 636 in Figure 1H); and when compared within roots (Figure 1I), or within January (Figure 1F), the number of unique OTUs in the roots of X792 in January was the highest (219 in Figure 1F; 189 in Figure 1I).

These data suggest that month, mulberry cultivar., and organ all contributed to the observed variation in the composition of the mulberry endophytic bacterial OTUs.

The obtained sequences were classified into 25 phyla, 56 classes, 169 orders, 325 families, and 647 genera. The bacterial composition and relative abundances varied across different samples. Figure 2 showed the diversity of bacterial communities in different samples at the phylum level. The dominant bacterial phyla were Proteobacteria and Actinobacteria across all samples; Actinobacteria was the predominant phylum in JSDS (45.88%), while Proteobacteria (47.65–88.44%) was the predominant phylum across all other samples. The relative abundance of Proteobacteria was higher in all cultivar X792 samples compared with the DS samples, except for the stem in February, while Actinobacteria showed the opposite trend. For Proteobacteria, Actinobacteria, and Firmicutes, there was no significant difference in relative abundance between mulberry cultivars; while there was a significant increase over relative abundance of Actinobacteria and Firmicutes in the stem of January compared with February (p < 0.01), and Proteobacteria showed the opposite trend, but not in root samples. Proteobacteria showed a significant increase in stems compared with roots in February (p < 0.01), and Actinobacteria showed the opposite trend, and Firmicutes showed significant increase in stems compared with roots in January (p < 0.01; Supplementary Table S3). These analyses indicated that the temperature had a greater influence on the endophytic bacterial content of the stem compared with that of the root.

Clustering of the top 30 genera was shown in Figure 3 (with additional supporting data present in Supplementary Table S4). These 30 bacterial genera belonged to five phyla including Proteobacteria (15 genera), Actinobacteria (12), Firmicutes (1), Bacteroidetes (1), and unclassified_k_norank_d_Bacteria (1). Pseudomonas, Steroidobacter, Rhodococcus, Ralstonia, Mycobacterium, and Cryptosporangium were the most abundant (>0.1%) genera across all samples, but the genera distributions differed greatly across different samples. Pseudomonas was the predominant genus in FSDS and FS792, while Steroidobacter was the most abundant genus in JRDS, JR792, FRDS, and FR792, and Rhodococcus was the predominant genus in JSDS and JS792. These results showed that the most predominant genus of the stem in January was Rhodococcus, and in February was Pseudomonas, while Steroidobacter was the predominant genus in the root, and there was no significant difference in the relative abundance of genus between mulberry cultivars. The clustering analysis showed a clear similarity among samples from each tissue (roots are similar and shoots are similar). For roots, two cultivars were gathered separately, and for stems, there were more similarities in January, and the same rule occurred in February.

The relative abundance of the top 15 core genera was also compared, and some differences were observed (Figure 4). The relative frequencies of Frigoribacterium (p = 0.04) and Pseudokineococcus (p = 0.04) were higher in JSDS compared with those in JS792 (Figure 4A), but there was no significant difference between JRDS and JR792 (Figure 4B), while that of the genus Gaiella (p = 0.013) was significantly lower in FSDS compared with that in FS792 (Figure 4C). The relative abundances of Pseudomonas (p = 0.004) and Acidibacter (p = 0.012) were significantly lower and significantly lower, respectively, in FRDS compared with those in FR792 (Figure 4D).

The relative abundance of Pseudomonas was also significantly lower in JSDS compared with that in FSDS (p = 0.0002), while that of Ralstonia was significantly higher in JSDS compared with that in FSDS (p = 0.005; Figure 4E). Moreover, the relative abundances of Pseudomonas (p = 0.000009) and Rhodococcus (p = 0.0007) were significantly lower and significantly higher, respectively, in JS792 compared with those in FS792, and that of Ralstonia was also significantly higher in JS792 compared with that in FS792 (p = 0.004; Figure 4F). Pseudomonas had a significantly lower abundance in JRDS compared with that in FRDS (p = 0.012; Figure 4G), and had a significantly lower relative abundance in JR792 compared with that in FR792 (p = 0.002; Figure 4H). The relative abundances of Rhodococcus (p = 0.003) and Ralstonia (p = 0.001) were significantly higher, and that of Steroidobacter was significantly lower, in JSDS compared with JRDS (p = 0.01; Figure 4I). In addition, the relative abundance of Rhodococcus was significantly higher in JS792 compared with that in JR792 (p = 0.0003), while that of Steroidobacter was significantly lower (p = 0.003), and those of Ralstonia (p = 0.002) and Pseudomonas (p = 0.008) were significantly higher, in JS792 compared with JR792 (Figure 4J). Furthermore, the genus Pseudomonas had a significantly higher relative abundance in FSDS compared with that in FRDS (p = 0.0002), while relative abundances of Steroidobacter (p = 0.006) and Virgisporangium were significantly lower in FSDS compared with those in FRDS (p = 0.006; Figure 4K). Finally, the relative abundances of Pseudomonas (p = 0.00002), Rhodococcus (p = 0.0005), and Ralstonia (p = 0.00005) were significantly higher, and that of Bradyrhizobium (p = 0.0009) was significantly lower, in FS792 compared with FR792, and that of Steroidobacter (p = 0.007) was significantly lower in FS792 compared with that in FR792 (Figure 4L). Overall, there were more significant differences in the relative abundances of the top 15 core genera between stem and root, especially in February, followed by the stem in January and February, while there were fewer significant differences between cultivars.

To further compare the relationship of endophytic bacteria populations among the stem and root of two mulberry cultivars in January and February, principal coordinates analysis (PCoA) based on Euclidean distances with arithmetic mean clustering was conducted using the genera. This analysis revealed the main variations in bacterial community composition and abundance among the samples. The PCoA results graphically demonstrated that organs or months were strong factors in accounting for the observed variations in the composition of the endophytic bacterial community, in which samples of the stem in January were placed at a higher PC 1 value (63.49%), while samples of the stem in February appeared a higher PC 2 value (27.54%), and root samples showed higher values of both PC 1 and PC 2 (Figure 5), and p = 0.001. Samples of the roots of different cultivars in different months clustered together, but the stems showed some separation, with samples of the stems in January and February clustering together, respectively. The PCoA results were supported by non-metric multidimensional scaling (NMDS) plots (p = 0.001; Supplementary Figure S2). In summary, these analyses revealed distinct differences in endophytic bacterial communities of stem and root, and stem of January and February, while no clustering was evident due to cultivars or root of January and February. This finding confirmed the results for the differential analysis of genera.

To explore the complexity of the interactions within the endophytic communities among the different samples, a correlation network analysis was conducted, and its topological properties were calculated. The correlation network analysis revealed that there was a difference between January and February. The complexity and modular structure were higher in January (Figure 6A) compared with that in February (Figure 6B), specifically, the average number of connections per node was higher in January (node average degree = 26.60) relative to the February samples (node average degree = 21.57; Table 2). Furthermore, January also presented a higher number of negative correlations (negative edges = 210) compared with February (negative edges = 113) on the basis of the same number of positive correlations. Additionally, the correlation network analysis of microbial communities for the sensitive or resistant cultivar also revealed differences. The sensitive mulberry cultivar DS (Figure 6C) appeared to have a more complex correlation compared with the resistant mulberry cultivar X792 (Figure 6D). The node average degree, positive edges, and negative edges in DS (node average degree = 22.07, positive edges = 163, and negative edges = 168) were all higher than those in X792 (node average degree = 15, positive edges = 120, and negative edges = 90).

The correlation network analysis was performed to assess the complexity of the interactions among the microbial taxa, the results also revealed a strong difference between the communities based on organs, such as stem (Figure 6E) and root (Figure 6F). The edge number for stem (edge = 218) was much higher compared with that of root (edge = 28), with 130 positive edges (positive edges of stem = 160, positive edges of root = 30) and 26 negative edges (negative edges of stem = 58, negative edges of root = 32). Simultaneously, the average number of connections per node was also much higher in the stem (node average degree = 15.03) compared with that in the root (node average degree = 4.43). In summary, the correlation network indicated that the interaction degree of the microbial community of mulberry was strongly influenced by the month, cultivar resistance, and organs. January, sensitive mulberry cultivar., and stem possessed a greater microbial complexity and abundance compared with February, the resistant mulberry cultivar., and root. The nodes with the highest connections in January were Quadrisphaera, Steroidobacter, Pseudokineococcus, Dokdonella, Streptomyces, Rhodomicrobium, Methylobacterium, Virgisporangium, Frigoribacterium, and unclassified_k__norank_d__Bacteria with a degree of 29, and in February was unclassified_p__Proteobacteria, unclassified_k__norank_d__Bacteria, and Pseudomonas with the degree of 26. Burkholderia-Caballeronia-Paraburkholderia, and unclassified_k__norank_d__Bacteria were the genera with the highest degree of 23 in cultivar X792, and Methylobacterium was the genus with the highest degree of 23 in cultivar DS. Ralstonia (degree = 10) and Pseudomonas (degree = 24) were the genera with the highest connections in root and stem, respectively.

Functions of microbial communities from all samples were predicted using PICRUSt2 on level 1 and level 2 (Figure 7; Supplementary Tables S10, S11). Genes associated with environmental information processing (level 1), membrane transport, and signal transduction (level 2) were significantly (95% CIs, p < 0.05) more abundant in the stem of resistant cultivar X792 compared with the sensitive cultivar DS in January, which may be related to low-temperature resistance, and warrants further research. Furthermore, Welch’s t-test results (Supplementary Table S6) indicated that, except for the environmental information processing of stems in January, there were no significant differences in several predicted pathways (e.g., metabolism, environmental information processing, cellular processes, genetic information processing, human diseases, and organismal systems) among the other samples based on cultivars. In addition, there were no significant differences between the root of X792 in January and February for the five pathways on level 1, as well as the JRDS and FRDS for genetic information processing and organismal systems pathways, but there were significant or significant differences in the other samples based on months. Except for the genetic information processing pathway of the JS792 and JR792 samples, there were no significant differences among the other samples in January based on organs, but there were significant or significant differences throughout the five pathways in all samples in February. The stems of January and February had different degrees of significant differences in the five pathways, as did the stem and root of February. Overall, Environmental information processing was the pathway exhibiting the most differences, followed by the cellular process.