PART 1 Revealing The Impact Of The Environment On Cistanche Salsa: From Global Ecological Regionalization To Soil Microbial Community Characteristics

ABSTRACT: To understand the regulatory relationship between the environment and Cistanche salsa, here we integrated the macro-and micro dimensions methods. From a macro perspective, the MaxEnt model indicated that countries along the Belt and Road Initiative, such as China, Egypt, and Libya, were particularly suitable for the growth of C. salsa from ancient times (Last Glacial Maximum and mid-Holocene) to the future (2050 and 2070). The Jackknife test revealed that precipitation is an important ecological factor that affects C. salsa’s distribution. From a macro perspective, 16S rRNA amplicon sequencing data showed that the soil microbial communities of three ecotypes (desert−steppe, grassland, and gravel−desert) were significantly different (p < 0.001). Core microbiome analysis demonstrated that the bacterial genera Arthrobacter, Sphingomonas, and Bacillus were enriched core taxa of C. salsa. LEfSe and random forest were used to excavate the Gillisia (desert−steppe), Flavisolibacter (grassland), and Variibacter (gravel−desert) as biomarkers that can distinguish among microbial communities from the three ecotypes. The prediction profile showed that the metabolic function of the microbial community was enriched in metabolic pathways and environmental information processing. Correlation analyses revealed that the altitude, precipitation of the warmest quarter (bio18), mean diurnal range (bio2), and mean temperature of the warmest quarter (bio10) were important ecological factors that affect the composition of soil microbial communities. This work provided new insights into the regulatory relationship among the suitable distribution of C. salsa, soil microbial communities, and ecological drivers. Moreover, it deepened the understanding of the interaction between desert plants and ecological factors in arid environments.

KEYWORDS: Cistanche salsa, MaxEnt, 16S rRNA amplicon sequencing, soil microbial community, environmental impact.

Cistanche deserticola has many effects, click here to know more

For more information please contact: Joanna.jia@wecistanche.com

1. INTRODUCTION

The economic value of medicinal plants has increased rapidly in recent years with the worldwide development and utilization of natural plants. The fleshy stem of Cistanche salsa is an edible and medicinal plant organ that is commonly used as a kidney tonic, an aphrodisiac, an antiaging and antioxidation treatment, an intestinal laxative, a liver protection treatment, and an anti-radiation treatment. Modern phytochemical studies on C. salsa have resulted in the identification and isolation of considerable quantities of bioactive compounds, such as various phenylethanoid glycosides, iridoids, alditols, and lignin, which have important medical and edible values.1 With the rapid development of the market for medicine and healthy foods, the demand for C. salsa has sharply increased and is accompanied by the overexploitation and plunder of wild resources. Therefore, the large-scale artificial introduction and cultivation of C. salsa have become an important measure to protect wild resources and slow ecological desertification. However, research on the artificial cultivation of C. salsa remains limited.

The environment has a massive impact on the growth and development of C. salsa.2 Plants grow in different environments, forming different ecotypes that exhibit significant changes inactive ingredients and gene expression.3,4 At a macro level, research has found that temperature, light, precipitation, and soil types all have an effect on plant growth and its active ingredients.5,6 The study of the interaction between the environment and plants and the comprehensive characterization of different ecotypes of plants are beneficial to the selection and cultivation of high-quality varieties.7 From a macro perspective, plants can have a considerable impact on soil organisms even if they have little or no direct contact with the soil system through their roots.8,9 However, current research on the relationship between plants and the environment, including climate factors and soil microbial community, is still unclear.

The species distribution model (SDM) is a statistical model that was established using existing species distribution data and environmental variables to infer the ecological needs of species and project their potential distribution areas,10,11 explore suitable growth areas by combining climate and soil factors, determine the appropriate environment for medicinal plants, and scientifically protect and grow endangered plants.12,13 SDMs, such as MaxEnt and biomod2, have successfully responded to the predictions of distribution trends imposed on endangered and ecological plants during climate change.14,15

16S rRNA amplicon sequencing of plant rhizosphere soil samples has been performed to explore the diversity of microbial communities, providing new insights into the relationship between plants and soil microbial communities.16,17

In this study, we explored the relationship between the environment and C. salsa from the macro-and micro dimensions. We performed the following. (1) We used the MaxEnt model to predict the global suitable growth areas of C. salsa from ancient times to the future (five periods: Last Glacial Maximum [LGM], mid-Holocene [MH], present, 2050, and 2070) and calculated the suitable areas of different levels and the contribution rate and range of bioclimatic variables affecting the distribution of C. salsa. (2) In combination with fieldwork, we collected three ecotypes (desert−steppe, grassland, and gravel−desert) soil samples from the best growing area (Tacheng and Xinjiang) of

C. salsa. We performed 16S rRNA amplicon sequencing to explore the characteristics of soil microbial communities. We also compared the differences in soil microbiomes in the three ecotypes and determined the core microbiomes and biomarkers that could distinguish among the three ecotypes. (3) We conducted correlation analysis and redundant analysis on the basis of the abundance of the core microbiomes, biomarkers, and bioclimatic variables to explore the regulatory relationship between C. salsa and the environment.

Cistanche has many effects in anti-fatigue

2. MATERIALS AND METHODS

2.1. Ecological Niche Modeling.

In this research, five resources were used to search for the present sites of C. salsa worldwide: (1) GBIF ; (3) National Specimen Information Infrastructure (4) published literature; and (5) fieldwork. Sample data without latitude and longitude were based on an online map latitude and longitude query The buffer analysis method was used to proofread and filter obtained distribution points to eliminate the effect of overfitting simulation caused by a large spatial correlation. The spatial resolution of the bioclimatic variables was 2.5 arc-min (about 4.5 km2), and the buffer distance was set to 3 km. Only one distribution point was retained when the distance between the distribution points was less than 3 km. After duplicate points were removed, 76 occurrences, all from the genuine production area of C. salsa, were collected and used for analysis (Figure 1a, Supplementary File 1).

We used the 19 bioclimatic variables of WorldClim as environmental predictors. The present bioclimatic variables data of this study were collected from the monitoring data of the WorldClim version 1.4 database18 during 1960−1990 with a resolution of 2.5 arc-min. The ancient (LGM and MH) and future (future 2050 and future 2070) bioclimatic variables available in WorldClim version 1.4 at 2.5 arc-min resolution were employed as predictors of the species distribution models. To avoid multicollinearity, we ran a correlation analysis on background points and eliminated one of the variables in each pair with a Pearson

correlation value > 0.819 (Figure S1). The eight bioclimatic variables finally included in the models were the mean diurnal range (bio2), the mean temperature of the wettest quarter (bio8), the mean temperature of the driest quarter (bio9), the mean temperature of the warmest quarter (bio10), annual precipitation (bio12), the precipitation of the driest quarter (bio17), the precipitation of the warmest quarter (bio18), and the precipitation of the coldest quarter (bio19).

We used the open-source software package maximum entropy model (MaxEnt v.3.4.0),20 which can be downloaded to build a species distribution model for C. salsa. The following parameters were used in the basic settings: random seed and a random test percentage of 25 and 10 replicates. By setting the random test percentage to 25%, we randomly selected 75% of the distribution points for the training set. By setting the number of replicates to 10, we ran the model 10 times with the same settings and averaged the output of all runs to obtain the final result. The area under the receiving operator curve (AUC) was used to evaluate the model’s goodness-of-fit, and the model with the highest AUC value was considered as the best performer. The jackknife procedure was used to assess the importance of the variables. Response curves were used to obtain the range of bioclimatic variables.

ArcGIS was used to analyze the bioclimatic variables that affect the distribution of C. salsa and to classify and calculate the area suitable for production.

2.2. Soil Sample Collection and Description. C. salsa

is naturally distributed in the three ecotypes of gravel−desert, grassland, and desert−steppe. In April 2017, we collected fleshy stem and soil samples representing the major ecotypes of C. salsa in Tacheng, Xinjiang, China (Figure 1b, Table 1). Gravel−desert samples were collected from Hejiaoke, Toli County (HJ1, HJ2, and HJ3). Grassland samples were collected from Yumin County (YM1, YM2, YM3, and YM4). Desert− steppe samples were collected from Jiang Alhan (JA1, JA2, JA3, JA4, JA5, and JA6). The soil samples we collected are all from the soil on the surface of C. salsa and its host parasitic site. Voucher specimens with

voucher numbers from 20170510079-DT to 20170510091-DT were deposited in the herbarium of the Institute of Medicinal Plant Development at the Chinese Academy of Medical Sciences in Beijing, China. After being cleaned, succulent stem tissues were cut into small pieces, immediately frozen in liquid nitrogen, and then stored at −80 °C until further processing. Soil cores were taken at a depth of 20 cm using

a stainless-steel cylindrical driller with a diameter of 5 cm and then stored at −20 °C in a portable refrigerator. After being transported to the laboratory, the soil samples were passed through a 2 mm sieve to remove plant tissues, roots, rocks, and other debris and then stored at −20 °C in a refrigerator before further experiments.

2.3. DNA Extraction and 16S rRNA Sequencing.

Soil DNA was extracted using a PowerSoil DNA Isolation Kit (MoBio Laboratories, Carlsbad, CA) in accordance with the manual. The purity and quality of the genomic DNA were checked on 0.8% agarose gels. The V3-4 hypervariable region of the bacterial 16S rRNA gene was amplified with the primers 338F (ACTCCTACGGGAGGCAGCAG) and 806R (GGACTACHVGGGTWTCTAAT).21 For each soil sample, a 10- digit barcode sequence was added to the 5′ end of the forward and reverse primers (Allwegene Co., Beijing). PCR was performed on a Mastercycler Gradient (Eppendorf, Germany) using 25 μL reaction volumes containing 12.5 μL of KAPA 2G Robust HotStart ReadyMix, 1 μL of forwarding primer (5 μM), 1 μL of reverse primer (5 μM), 5 μL of DNA (total template quantity is 30 ng), and 5.5 μL of H2O. The cycling parameters were as follows: 95 °C for 5 min, followed by 28 cycles of 95°C for 45 s, 55 °C for 50 s, and 72 °C for 45 s, with a final extension at 72°C for 10 min. Three PCR products per sample were pooled to mitigate reaction-level PCR biases. The PCR products were purified using a QIAquick Gel Extraction Kit (QIAGEN, Germany) and then quantified using real-time PCR. Deep sequencing was performed on a MiSeq platform by Allwegene Co. (Beijing). After the run, image analysis, base calling, and error estimation were performed using Illumina Analysis Pipeline Version 2.6.

2.4. Analysis of 16S rRNA Amplicon Sequencing Data.

All sequencing data were submitted to the NCBI Short Archive (SRA) under the SRA submission SUB7456002. Raw data were first screened, and sequences were removed on the basis of the following considerations: sequences shorter than 200 bp with low-quality score (≤20) and containing ambiguous bases or not matching primer sequences and barcode tags. Qualified reads were separated using sample-specific barcode sequences and trimmed with Illumina Analysis Pipeline Version 2.6. Then the data sets were analyzed using QIIME. The sequences were clustered into operational taxonomic units (OTUs) at a similarity level of 97%22 to generate rarefaction curves and calculate richness and diversity indices. The Ribosomal Database Project Classifier tool was used to classify all sequences into different taxonomic groups, in which the confidence threshold is set to 0.7.23 Clustering analyses were performed on the basis of OTU information from each sample by using R 3.6.1 to examine the similarity between different samples.24 The UniFrac distances matrix between microbial communities from each sample were calculated using the Tayc coefficient and represented as an unweighted pair-group method with arithmetic mean clustering tree, which describes the dissimilarity (1- similarity) among multiple samples.25 A Newick-format tree file was also generated through this analysis. Alpha diversity was applied to analyze the complexity of species diversity for a sample using four indexes, namely, Chao1, observed species, and Shannon and Fisher

diversity indexes. These indexes were calculated using QIIME software (Boulder, CO, USA) in Python (v.1.8.0) (La Jolla, CA, USA).26 Beta diversity analysis was used to evaluate differences in samples in terms of species complexity. Beta diversity was calculated using the principal coordinate analysis (PCoA) and cluster analysis in QIIME.27 Analysis of molecular variance (AMOVA) was performed using mother.28 The Kruskal−Wallis test was used to calculate the OTU difference between groups (p-value <0.05), and heat maps were drawn using pheatmap (R package). Core microbiome analysis was adopted from the core function in the R package microbiome (sample prevalence = 20%, relative abundance = 0.01%) by MicrobiomeAnalyst29 (https://www. microbiomeanalyst.ca/MicrobiomeAnalyst/home.xhtml).

Linear discriminant analysis (LDA) and random forest (RF) methods in the MicrobiomeAnalyst Web site were used to determine microbiome biomarkers. A nonparametric factorial Kruskal−Wallis sum−rank test was first performed to identify features with significant differential abundance considering the experimental factor or class of interest. Next, LDA (the threshold is set to 2) was performed to calculate the effect size of each differentially abundant feature. Features were considered significant on the basis of their adjusted p-value. The default adjusted p-value cutoff was 0.05. RF analysis was performed using randomForest package5. This method utilizes an ensemble of classification trees, each of which is grown via random feature selection from a bootstrap sample at each branch.

Tax4Fun (R package, http://tax4fun.gobics.de/) was used to predict the microbial functional profiles of microbiomes in the soil samples. The OTU Biom table of the soil microbiome was used as an input file for the metagenome imputation of C. salsa soil samples. Then the predicted gene class abundances were analyzed at the KEGG Orthology (KO) group level 3. The results from Tax4Fun were analyzed in Doby (R package).

2.5.Correlation Analysis of Key Microbial Communities and Bioclimatic Factors.

We used ArcGIS to numerically extract key bioclimatic factors from the 13 soil sampling points. Redundancy analysis of key genera of the microbial communities (five biomarkers, five core microbiomes) and bioclimatic factors was performed using

Canoco 5 software. Log2 data conversion was uniformly performed before the analysis. Spearman correlation coefficients were calculated for the abundance of five biomarkers, five core microbiomes, and bioclimatic factor data integration by applying SPSS. The correlation analysis results were drawn by pheatmap (R package).

Cistanche can resist aging

3. RESULTS

Prediction of the Global Distribution Shift of C. salsa in the Different Periods. The calculated ROC showed

that the AUC value was 0.977 (Figure S2), indicating that the model agreed well with the data.30 Habitat suitability simulation results for five periods (LGM, MH, present, 2050, and 2070) were illustrated in Figure 2 and Table 2. C. salsa-growing suitability was divided into five categories on the basis of statistical thresholding: unsuitable (class 5: 0−20%), marginal (class 4: 20−50%), fair (class 3: 50−75%), good (class 2: 75−

90%), and excellent (class 1: 90−100%). This article only discussed the suitable distribution areas for class 1, class 2, and class 3. The present potential distributions of C. salsa indicated that this species was widely distributed across three climatic zones: temperate desert climate, subtropical Mediterranean climate, and tropical desert climate regions. All continents, except Antarctica, contain regions that were suitable for C. salsa. These regions were mainly distributed in central and western Asia and northern Africa and were scattered in central and western North America, central South America, and western Oceania. Excellent (class 1) suitable areas were the most widely distributed in Egypt (76 750 km2), whereas the good (class 2) and fair (class 3) suitable areas were most widely distributed in China (class 2: 722 075 km2 and class 3: 1 024 600 km2). The ecologically suitable areas in Asia were primarily restricted to China (class 1: 45 775 km2), Jordan (class 1: 16 075 km2), Israel (class 1: 14 975 km2), Saudi Arabia (class 1: 14 925 km2), and Iran (class 1: 12 600 km2). The suitable areas of C. salsa in Africa were mainly distributed in Egypt, Libya (class 1: 34 400 km2), and Tunisia (class 1: 275 km2). Suitable areas for C. salsa in South America were mainly distributed in Chile (class 1: 16 550 km2).

From a spatial perspective (Figure 2), the suitable regional change trend for the five periods first increased and then decreased with the largest area in 2050 (Figure 3a, Figure S3). We saw a large decrease in suitable C. salsa area within the LGM relative to that in the other four periods. Compared with the present (4 358 775 km2), suitable areas decreased by 50% (2 160 975 km2) in the LGM period, by 10% (MH: 3 910 350 km2) in the MH period, and by 1% (4 328 800 km2) in the 2070 period and increased by 2% (4 428 950 km2) in the 2050 period. Notably, from the present to the future (from 1960 to 2080), the excellent area (class 1) for C. salsa gradually shrank (present, 243 200 km2; 2050, 232 425 km2; 2070, 110 800 km2).

Cistanche extract has many health benefits

3.2. Prediction of the Suitable Distribution Areas of C. salsa in China and the Range of Bioclimatic Variables.

The suitable areas for C. salsa growth in China were mainly restricted to northwestern Xinjiang, central Inner Mongolia, northern Shaanxi, northern Shanxi, northeast Qinghai, northern Gansu, and most of Ningxia (Figure 2b). The most suitable areas in China were mainly distributed in Xinjiang, of which the excellent area (class 1) covered 23 650 km2, the good area (class

2) covered 220 950 km2, and the fair area (class 3) covered 189 975 km2. Among the good areas (class 2), Xinjiang accounted for 60.99%, followed by Inner Mongolia (86 300 km2) and Gansu (24 450 km2) with 23.82% and 6.75%, respectively (Figure 3b). The results of the MaxEnt model indicated that Xinjiang was particularly suitable for C. salsa growth. However, different habitats, such as gravel−desert, grassland, and desert−steppe, were found in the same area during field trips. Therefore, we were highly interested in the characteristics of soil microbial communities and their relationship with the environment in different ecotypes of C. salsa in Tacheng, Xinjiang, China.

The importance of each bioclimatic variable to the distribution of C. salsa is shown in Table 2. All values were means of 10 replicate runs. The top three bioclimatic variables with the greatest influence on the C. salsa distribution were the precipitation of the driest quarter (bio17), the precipitation of the warmest quarter (bio18), and the mean diurnal range (bio2) with relative contribution rates of 25.81%, 17.65%, and 13.70%, respectively. As indicated by the response curve obtained by the MaxEnt model, the range of the bio factor could be calculated when the prediction probability exceeded 0.5. For example, the range of bio17 was from −135.30 to 11.34 mm, that of bio 18 was from −206.60 to 4.40 mm, and that of bio2 was from 11.65 to 12.50 °C.

3.3. Soil Microbiomes From Three Ecotypes of C. salsa Exhibiting Distinct and Overlapping Microbial Communities.

16S rRNA sequencing resulted in 518 217 raw reads, among which 441 576 were screened for quality and length. The data set contained 11 818−26 431 (mean: 20 150) sequences per sample. The high-quality reads were clustered on the basis of

>97% sequence identity into 2 788 microbial OTUs (Table 1).

The microbial community was classified into 34 phyla and 321 genera. At the phylum level (Figure 4a), Actinobacteria (DS, 31.94%; GD; 46.42%; GL, 33.33%), Proteobacteria (DS,

23.25%; GD, 22.53%; GL, 24.68%), and Gemmatimonadetes (DS, 17.77%; GD, 8.02%; GL, 8.36%) were dominant in the three ecotypes. At the genus level (Figure 4b), the desert− steppe was dominated by Euzebya (4.82%) and Arthrobacter (1.74%), whereas the most dominant genera in the gravel− the desert was Arthrobacter (8.35%) and Bacillus (4.95%). Bacillus (6.89%) and Arthrobacter (5.57%) were dominant in the grassland. The top 10 abundant microbial communities in the desert−steppe were classified into eight phyla (Figure 4c), those in the gravel−desert were classified into seven phyla (Figure 4d), and those in the grassland were classified into 10 phyla (Figure 4e).

The measurements of within-sample diversity (α-diversity) revealed a diversity change from the grassland to the gravel− desert and desert−steppe (Figure 4f). The α-diversity of soil microbial communities in each sample was evaluated on the Shannon, Chao 1, Fisher, and observed species diversity indexes. The Shannon, Chao 1, observed species, and Fisher indexes suggested that the α-diversity of the grassland soil communities

was higher than that of the other two ecotype communities. The results of the rarefaction curves (Figure S5) were similar to the above results, except for YM1 and HJ3. AMOVA results (Table S3) showed that there are significant differences (p < 0.01) among the three ecotypes. The results of the unconstrained PCoAs of unweighted UniFrac distance 2D plots (Figure 4g) indicated that the soil samples of bacterial microbiomes from different ecotypes, except YM1 and HJ3, were well clustered. The Bray distance diversity tree clustering based on single-clustering algorithm results (Figure 4h) of the three ecotypes soil samples demonstrated that the grassland samples, except YM2, were closely clustered, the desert−steppe soil samples, except JA1, were closely clustered, and the gravel−desert soil samples were closely clustered.

3.4. Differential Microbiome Screening of Three Ecotypes of C. salsa.

The heat maps of the abundances of different genera (Figure 5a) indicated that the genus abundance of samples from the desert steppe was different from that of samples from the other two ecotypes. As shown in the results of

LEfSe (Figure 5c) and RF (Figure 5d), the mini heat map to the right indicates the abundance of the microbial features in the three ecotypes at the genus level. The genera that can represent the desert−steppe ecotype include Gillisia, Illumatobacter, Salegentibacter, Marinimicrobium, etc. Variibacter is a symbolic

genus level that can represent the gravel−desert ecotype. The biomarkers in the grassland ecotype contain Flavisolibacter and Agromyces. After the results of the two methods were combined, 11 biomarkers were selected (Table S4). Figure 5b shows the heat map of the abundance of 11 biomarkers.

3.5. Core Microbiome Screening and Metabolic Function Prediction of the Three Ecotypes of C. salsa.

A Venn diagram (Figure 5e) was plotted for the OTUs obtained from all soil samples, and the results showed that the three ecotypes shared 1712 OTUs. The persistence method was adopted from the core function in the R package microbiome to identify the core microbiome in the three ecotypes of C. salsa. This core bacterial microbiome contained six OTUs and corresponded to 19.64% of the whole microbiome. With the exclusion of undefined and duplicated genera, these OTUs were classified into six genera, and their abundances were drawn in a heat map (Figure 5f).

The functional prediction results (Figure 6, Supplementary File 2) suggested that the functional metabolism pathways of soil microbiomes in the three ecotypes of C. salsa were identical in carbohydrate and that amino acid metabolisms were abundant among the metabolic pathways. Membrane transport and signal transduction were also abundant in environmental information processing.

3.6.Correlation Analysis between Microbial Communities and Bioclimatic Variables of Three Ecotypes of C. salsa.

The redundant analysis of the core, biomarker microbiome abundance, and bioclimatic variables was performed at the genus level, and reanalysis was performed on the basis of effects. The adjusted interpretation of variance was 32.5%. Precipitation of the warmest quarter (bio18) explained 23.9% of microbial communities (p = 0.07). Mean temperature of warmest quarter (bio10) and mean diurnal range (bio2) were positively correlated with Illumatobacter and negatively correlated with Bacillus (Figure 7a).

Correlation analysis was conducted for the core, biomarker microbiome abundance, and seven bioclimatic variables. The correlation network results (Figure 7b) revealed that Illumatobacter and Salegentibacter (biomarkers in the desert−steppe) were significantly positively correlated with the mean diurnal range (bio2) and mean temperature of the warmest quarter (bio10) but significantly negatively correlated with altitude (alt) and precipitation of the warmest quarter (bio18). On the contrary, Agromyces (biomarker in the grassland) was positively correlated with altitude (alt) and precipitation of the warmest quarter (bio18) but was negatively correlated with mean diurnal range (bio2) and mean temperature of the warmest quarter (bio10). In addition, Arthrobacter (core microbiome) was significantly negatively correlated with the mean diurnal range (bio2) and the mean temperature of the warmest quarter (bio10). Rubrobacter was significantly positively correlated with annual precipitation (bio12).