Evaluating the Endangered Mechanism and Different Conservation Effects of Pinus squamata from the Rhizosphere Soil Bacterial Community Perspective

1. Introduction

China is home to a high diversity of gymnosperms [1], and the Chinese flora includes many ancient relict gymnosperms species, including Ginkgo biloba, Metasequoia glyptostroboides, Glyptostrobus pensilis, Cathaya argyrophylla, Pseudolarix amabilis and others. However, many species are severely threatened, with habitat degradation, restricted distribution, and over-exploitation listed as the top three threats to endangered gymnosperm species [2]. Pinus Sect. Strobus is thought to have evolved in China [3]. The section is distinguished by the needles, which are in clusters of five, and by the single-vascular bundle. Sect. Strobus is therefore considered to be the earliest-diverged and most ancient lineage of the Pinaceae. The threats facing the species of Pinus Sect. Strobus and the circumstances underlying their scarcity have been studied [4].

Pinus squamata, also known as the Qiaojia five-needle pine, was published as new species in 1992 [5]. The tree has mottled bark and a straight trunk, and has high ornamental and economic value. Currently, only 34 individuals of P. squamata are known from the wild. The species has been assessed as Critically Endangered (CR) according to the endangered species classification system adopted by the International Union for Conservation of Nature and Natural Resources (IUCN) in 2001. Indeed, at the 24th World Congress on Conservation, the IUCN included P. squamata on the list of the 100 most endangered species in the world [6].

The Plant Species with Extremely Small Populations (PSESP) are defined as having small remaining populations, a restricted narrow habitat, severe anthropogenic disturbance and being at a high risk of extinction [7,8]. It should be emphasized that naturally rare species do not qualify as PSESPs, and external disturbance is necessary to determine whether a species is included as a PSESP [9]. Several national and provincial level programs and projects aiming to rescue PSESPs and supporting their conservation have been implemented [10,11], and the PSESP concept and the conservation action plans in China have attracted widespread attention in the field of conservation biology [12,13,14]. P. squamata was listed as one of 20 key species of PSESP in Yunnan and was also included on the Chinese national list of 120 PSESPs in 2011 [7,15].

The rhizosphere, which is the soil directly surrounding plant roots, is an important zone for the interactions between plants, soil and microorganisms [16]. Rhizosphere soils provide natural microhabitats for diverse microorganisms and are therefore considered to be hotspots of microbial diversity and activity in soils [17]. The rhizosphere microbiome affects plant growth and development, as well as resistances to stress and efficiency of nutrient use, therefore playing an important role in promoting soil quality and plant health [18,19]. Plant growth-promoting rhizobacteria (PGPR) are soil bacteria that are able to colonize the rhizosphere and to improve plant growth and health [20]. It have a crucial role in soil quality and fertility, as well as in the management of abiotic and biotic stresses [21]. PGPR can interact with plants directly, for example, by increasing the availability of essential nutrients, and can also indirectly affect plants by protecting them against diseases via competition with pathogens [22]. PGPR are under consideration as biostimulants for sustainable agriculture, alleviators of abiotic stresses in soil, and green bioinoculants and biofertilizers [23,24,25].

With the drastic changes occurring in the environment, soil biodiversity is variational in a rapidly changing world [26]. Detrimental condition of underground may cause species phenotypic decline, which may be causes for endangered. Therefore, the integrated study of soil properties and rhizosphere bacteria can assist conservation, as it allows us to understand the conditions required for conservation of rare and endangered species. Comparative analysis of rhizospheric fungi between wild, ex situ, and reintroduced P. squamata was researched by [27]. However, until now, the bacterial communities associated with the rhizospheric soil of P. squamata have remained unknown.

We hypothesized that the endangered cause of wild P. squamata populations is lower bacterial diversity in the rhizosphere soils, fewer beneficial taxa, and low soil nutrients. Ex-situ and reintroduced populations have higher bacterial diversity, beneficial taxa. Different populations vary widely in bacterialcomposition, so rhizosphere bacteria can be used to assess the effect of conservation and guide conservation practices. It is our hope that this study will inform the further optimization of conservationapproaches and strategies, and relevant plant growth-promoting bacteria (PGPR) will guide research intomanaging, protecting and restoring P. squamata habitat as well as further species conservation.

2. Materials and Methods

2.1. Sample collection

Samples were gathered from the wild (W), ex-situ (Ex) and re-introduced (Re) populations at different conservation sites in June 2020. The wild collection site was at Qiaojia County, Zhaotong City (WQ, 26°52′03.96″ N, 103°00′42.44″ E, 2206 m), ex-situ samples were collected from Caojian forestry farm, Dali Bai Autonomous Prefecture (EC, 25°45′32.73″ N, 99°06′53.94″ E, 2502 m), Yipinglang forestry station, Chuxiong Yi Autonomous Prefecture (EY, 25°08′08″ N, 101°54′01″ E, 1893 m), Kunming Botanical Garden (EK, 25°08′40.13″ N, 102°44′28.96″ E, 1990 m) and Qiaojia County (EQ, 27°0' 30.34" N, 102°57' 25.47" E, 1876.62 m). The samples from re-introduced population were taken from a site in Qiaojia County (RQ, 26°52′00″ N, 103°00′39.9″ E, 2133 m). Samples were collected from three individuals at each site. A total of 18 rhizospheric and 18 bulk soil samples were collected, respectively.

To collect a sample, we first removed the surface litter with a small shovel, gently dragging the shovel along the root to the tip, at a depth of 0 cm - 20 cm from the surface. Next, the soil only loosely bound to the root was shaken off, and was collected and labeled as bulk soil. Then, the soil within 1 mm - 10 mm of the root surface in four directions was collected into a sterile ziplock bag as rhizosphere soil, and was stored in liquid nitrogen [28]. Immediately on returning to the laboratory, the samples were examined: if the soil was too sticky, we used tweezers cooled on ice to remove impurities such as twigs, fine roots, and crushed stones. Samples were then mixed well, and those that could be sieved were passed through a 0.355 mm mesh sieve. Samples were then put into 5 mL cryopreservation tubes and stored in a -80 °C freezer.

2.2. Soil physical and chemical properties

Bulk soils were used to determine soil physicochemical properties. National standards of the People’s Republic of China were used to determine the soil physicochemical parameters [29]: soil pH (NY/T1377-2007), soil organic matter (OM, NY/T1121.6-2006), soil total nitrogen (TN, NY/T1121.24-2012), soil total phosphorus (TP, NY/T88-1988), soil total potassium (TK, NY/T87-1988), soil available nitrogen (AN), soil available potassium (AK, NY/T889-2004), soil available phosphorus (AP, NY/T1121.7-2014). Soil physical and chemical properties were tested by the Yunnan Sanbiao Agriculture and Forestry Technology Co., Ltd.

2.3. DNA extraction, PCR amplification, and high -throughput sequencing

Total DNA from the soil microbes was directly extracted using the Power Soil DNA Isolation Kit (MoBio Laboratories, San Diego, CA, USA). DNA purity and concentration were assessed with NanoDrop 2000 spectrometer (Thermo Fisher Scientifc, Wilmington, DE, USA), and DNA integrity was assessed using 1% agarose gel electrophoresis. The bacterial universal primer pairs 338F (ACTCCTACGGGAGGCAGCAG)_806R (GGACTACHVGGGTWTCTAAT) were then used to amplify the V3-V4 region of the 16S rRNA genes [30].

PCR amplification was conducted using TransGen AP221-02: TransStart FastPfu DNA Polymerase (TransGen Biotech, Beijing, China) and was performed on a GeneAmp 9700 thermal cycler (Applied Biosystems, Foster City, CA, USA). The reaction mixture included 4 µL of 5×FastPfu buffer, 2 µL of 2.5 mM dNTPs, 0.4 µL of FastPfu polymerase, 0.8 µL of the forward primer (5 µM), 0.8 µL of the reverse primer (5 µM), 0.2 µL of BSA, and 10 ng of template DNA, and double distilled water (ddH2O) was added to a final volume of 20 µL. Thermal cycling conditions were as follows: 3 minutes denaturing step at 95 °C, then 27 cycles (30 seconds at 95 °C; 30 seconds for annealing at 55 °C; and 45 seconds at 72 °C), 10 minutes for a final extension at 72 °C, 10 °C until halted. PCR amplification was detected using 2 % agarose gel electrophoresis, and the target fragments were cut out and recovered.

Next, the products were purified, quantified and homogenized to form a sequencing library, and the constructed library was checked for quality. Finally, qualifying products were subjected to bidirectional high-throughput sequencing using Illumina MiSeq PE300 (Illumina, San Diego, USA). The sequencing of all samples in this study was performed by the Shanghai Majorbio Bio-pharm Technology Co., Ltd. The National Microbiology Data Center (NMDC, https://nmdc.cn/) allows a huge amount of microbiological data to be organized and integrated in an effective way, and shared in an open manner [31]. The raw sequences of 18 samples have been deposited in NMDC with BioProject ID NMDC10018203 and accession numbers NMDC40026266- NMDC40026283.

2.4. Data analysis

FLASH (version 1.2.11) software was used to join the reads from each sample, and to obtain high-quality clean reads [32]. The software QIIME (version 1.9.1) and Fastp (version 0.19.6) were used to filter the clean tags and obtain effective tags. The sequences were then clustered into operational taxonomic units (OTUs) with UPARSE (version 7.0.1090) based on a 97% nucleotide similarity threshold [33]. Taxonomic assignments were performed using the RDP classifer algorithm and the Silva 16S rRNA database (version SSU 138), with a confidence threshold for taxonomic assignment set to 70% [34]. OTUs identified belonging to chloroplasts and mitochondria were removed from the dataset.

Alpha diversity metrics reflecting the richness and diversity of the communities were calculated using Mothur (version v.1.35.1). Beta diversity was also estimated as a representation of the compositional differences between communities. Principal coordinates analysis (PCoA) using the Bray-Curtis distance metric was used to evaluate similarities across community structures. Analysis of similarities (ANOSIM) was used to test whether the differences between groups was significantly greater than that within groups. Common and unique taxonomic communities among the different groups were visualized using a Venn diagram. Significant differences in relative abundance were tested with Kruskal-Wallis rank sum tests. Linear discriminant analysis Effect Size (LEfSe) measurements analysis was conducted to identify taxonomic biomarkers for different groups. The relationships between soil physicochemical parameters and bacterial communities were analyzed using redundancy analysis (RDA) and correlation heatmap analysis. Potential functions were predicted using the Phylogenetic Investigation of Communities by Reconstruction of Unobserved States (PICRUSt) program [35]. Functional Annotation of Prokaryotic Taxa (FAPROTAX), which is a tool for predicting ecologically relevant functions of bacterial and archaeal taxa derived from 16S rRNA amplicon sequencing [36]. Majorbio Cloud is a one-stop, comprehensive bioinformatic platform for multi-omics analyses [37], and our bioinformatics analyses were performed using the online Majorbio Cloud Platform (www.majorbio.com).

Statistical analysis was conducted using SPSS 22.0 software (SPSS Inc., Chicago, IL, USA) and results are given as means ± SD (standard deviations), with a P < 0.05 considered to be statistically significant.

3. Results

3.1. Sequencing quality

18 samples were subjected to sequence whole community 16S rRNA gene sequencing on an Illumina MiSeq PE 300 platform. After filtering out low-quality reads, we obtained 969,613 clean reads, with an average sequence length of 414 bp (Supplementary Table S1). With increasing numbers of sequencing reads, the rarefaction curves eventually became flat (Supplementary Figure S1), indicating that the sequencing depth of all samples is reasonable, and can truly and comprehensively reflect the structure and composition of the P. squamata rhizosphere bacterial community at different sites. In addition, the coverage of bacteria of samples was higher than 97.00% (Table 1), which also indicates that the integrity of the sequencing data is reliable, and that the probability of there being undetected microbial sequences in each sample was extremely low.

3.2. Alpha diversity

The alpha diversity at the OTU level is presented (Table 1). The sobs, ace, and chao1 indices were all highest at the RQ site, with values of 2404.67, 2999.92, and 2996.75, respectively. The EY site had the lowest values of these three indices, and these values were significantly different at EY from those of the other five sites (P < 0.05). The EY site thus had the lowest community richness and the RQ site had the highest. The Shannon index at the RQ site was the largest, with a value of 6.46, which was significantly different from that at either the EY or WQ sites (P < 0.05). The Simpson indices of the EK, RQ and EQ sites were smaller, with values of 0.0041, 0.0050, and 0.0059 respectively, and were significantly different from those at the EC or WQ sites (P < 0.05). These data suggest that the RQ site had high community diversity, while the WQ site had low community diversity.

3.3. The bacterial community in P. squamata rhizosphere soil

The bacterial communities of the 18 tested P. squamata rhizosphere soils comprised a total of 37 phyla, 119 classes, 281 orders, 431 families, 816 genera, 1786 species, and 5945 OTUs (Figure 1). Phyla with a relative abundance of less than 0.01 (1%) in all samples were classified as “others”. 13 phyla varied among the six sites (Figure 1A). The relative abundances of Actinobacteriota (20.42 % - 34.51 %), Proteobacteria (21.51 % - 39.42 %), Acidobacteriota (10.43 % - 23.41 %), Chloroflexi (5.37 % - 21.36 %) and Firmicutes (1.95 % - 4.59 %) accounted for 83.57 % (RQ) - 91.86 % (EY) of the total abundance. The relative abundance of the Actinobacteriota, Proteobacteria and Acidobacteriota exceeded 10 % of the total, and these phyla make up the dominant bacterial community in the P. squamata rhizosphere soil at all sites.

At the genus level, the composition of the bacterial community also differed among the P. squamata rhizosphere soil samples. Genera with a relative abundance of less than 0.05 (5 %) in all samples were categorized as “others” (Figure 1B). Bradyrhizobium (7.50 %) was abundant at the EC site. Acidothermus (9.85 %), Bradyrhizobium (5.18 %) and Mycobacterium (5.12 %) were dominant at the EY site. Bradyrhizobium (7.23 %) and Mycobacterium (5.91 %) were abundant at the WQ site. The relative abundance of Bradyrhizobium at EK, RQ, and EQ was 1.03 %, 3.64 % and 1.83 %, respectively.

Common and unique communities at the different sites were visualized with a Venn diagram (Figure 2). At the phylum level, 20 (54.05 %) of the 37 phyla were common to all of the six sites and were therefore defined as core phyla (Figure 2A). Soil from the EC site harbored 2 unique phyla. At the genus level, 193 (23.65 %) of the 816 genera were common to all the different sites, and were therefore considered to be core genera (Figure 2B). The number of genera unique to each site were EY (3), EQ (5), WQ (10), RQ (25), EK (30) and EC (36).

3.4. Beta diversity

Principal co-ordinates analysis (PCoA) of the bacterial communities of the P. squamata rhizosphere at the OTU level were presented (Figure 3). Samples were scattered (R = 0.93, P = 0.001), indicating that there were significant difference between the different sites. PC1 explained 43.53 % of the variance, and PC2 accounted for 16.83 %, together accounting for 60.36 % of the difference in bacterial community structure among sites.

3.5. Difference analysis of bacterial communities in P. squamata rhizosphere soil

Significant differences between relative abundances of the top 10 phyla and genera were tested with Kruskal-Wallis rank sum tests (Figure 4). At the phylum level, Gemmatimonadota was present at significantly different abundances at different sites (0.001 < P ≤ 0.01), while the abundances of Actinobacteriota, Proteobacteria, Chloroflexi, Myxococcota, Methylomirabilota, and Bacteroidota were different between sites (0.01 < P ≤ 0.05) (Figure 4A). Furthermore, at the genus level, Mycobacterium abundance differed between sites (0.01 < P ≤ 0.05), Bradyrhizobium and Acidothermus abundances were significantly different among different sites (0.001 < P ≤ 0.01) (Figure 4B).

LEfSe analysis was performed from phylum to genus level based on the all-against-all strategy. Only taxa above the linear discriminant analysis (LDA) significance threshold > 4.0 were presented in the six soil groups. At the phylum level, we found that the Proteobacteria were enriched at EC, the Chloroflexi at EY, the Bacteroidota at EK, and the Gemmatimonadota at RQ. The Myxococcota and Verrucomicrobiota were enriched at WQ, while the Actinobacteriota and Methylomirabilota at EQ. At the genus level, Bradyrhizobium was enriched at EC, Acidothermus and Conexibacter at EY, and Mycobacterium and Streptomyces at WQ (Supplementary Figure S2).

3.6. Relationships between bacterial communities and soil physicochemical factors

Redundancy analysis (RDA) was conducted to evaluate the relative effects of soil physicochemical factors on variation in the bacterial communities in P. squamata rhizosphere soil (Figure 5). At the phylum level, the X-axis and Y-axis were explained 29.80% and 21.38% of the variation in bacterial community composition, respectively (Figure 5A). Of the eight soil physicochemical parameters, soil TK and AN affected bacterial community structure (0.01 < P ≤ 0.05), while the other parameters did not (P > 0.05) (Supplementary Table S2). At the genus level, RDA1 explained 51.38%, and RDA2 interpreted 13.60%, together accounting for 64.98% of the total variation in bacterial community structure between sites (Figure 5B). pH and TP were the main factors extremely significantly influencing the composition of the bacterial community (P≤0.001). TK significantly affected community structure (0.001 < P ≤ 0.01), while AP and AK influenced community structure (0.01 < P ≤ 0.05). However, OM, TN, and AN did not affect community structure (P > 0.05) (Supplementary Table S3).

Correlations between soil physicochemical parameters and the relative abundances of the top 10 most abundant bacterial phyla were visualized using heatmaps (Figure 6). TK and other factors clustered into two branches. TK was significantly positively correlated with the Proteobacteria, and negatively correlated with the Methylomirabilota. TP was extremely significantly positively correlated with the Gemmatimonadota, and positively correlated with the Methylomirabilota. pH was significantly positively correlated with the Myxococcota and Methylomirabilota, positively correlated with the Bacteroidota, and negatively correlated with the Planctomycetota. AP was significantly positively correlated with the Methylomirabilota and Bacteroidota, and positively correlated with the Gemmatimonadota and Myxococcota. AK was positively correlated with the Actinobacteriota. OM was positively correlated with the Bacteroidota. AN and TN were negatively correlated with the Planctomycetota.

3.7. Function prediction in rhizosphere bacteria of P. squamata

PICRUSt2 was used to predict the functions of the bacteria in the P. squamata rhizosphere soil. COG function classification of the organisms in the P. squamata rhizosphere was mainly related to “function unknown”, “amino acid transport and metabolism”, and “energy production and conversion” (Supplementary Figure S3). The functions of the top 20 P. squamata rhizosphere bacteria (in total abundance) were analyzed using FAPROTAX, and a functional heatmap was obtained (Supplementary Figure S4). “Chemoheterotrophy”, “aerobic_chemoheterotrophy”, “nitrogen fixation”, “cellulolysis”, “nitrate_reduction” and “other functions” were found to be most abundant. Functional difference analysis between groups at different sites showed that there were significant differences in “cellulolysis” and “nitrogen fixation” (0.001 < P ≤ 0.01), no difference in “ureolysis” (P > 0.05), and differences in seven other types (0.01 < P ≤ 0.05) between sites (Supplementary Figure S5).

4. Discussion

4.1. Diversity of P. squamata rhizosphere soil bacteria

Plant-microbiome interactions are directly related both to microbial community assembly and to plant health [38]. In our study, the RQ site had the highest bacterial community richness and diversity, this site may therefore be suitable for furture reintroduction of P. squamata. Previous study indicated that the EK had the highest fungal community richness and diversity [27]. Of the four ex-situ sites, EK showed the highest rhizosphere bacterial diversity and was similar to EQ, indicating that the ex-situ P. squamata conservation projects taking place at EK could be expanded. Lower bacterial community richness was observed at EY, and the WQ and EY sites had lower observed bacterial diversity, which may be related to the older ages of the trees sampled at these sites.

4.2. Structure of the P. squamata rhizosphere bacterial community

In our study, 13 bacterial phyla varied among the different sites. The Actinobacteria, Proteobacteria, Acidobacteria and Chloroflexi were the dominant phyla in the bacterial communities of the P. squamata rhizosphere soil. A total of 12 bacterial phyla were annotated in the rhizosphere of P. dabeshanensis, of which the dominant phyla were the Proteobacteria, Acidobacteriota, Actinobacteriota and Chloroflexi [39]. However, a further phylum, the Patescibacteria, was unique to P. squamata. Patescibacteria have also been isolated from the rhizosphere of the salt-tolerant Suaeda salsa (Amaranthaceae) [40], and the presence of these bacteria may be related to plant salt or drought tolerance.

Actinobacteria, which can produce many biologically active compounds and degrade cellulose [41], have been marketed as being able to promote soil and plant health [42]. The relative abundance of Actinobacteria at Qiaojia was higher in soil from trees in ex-situ (EQ) conservation than soils from the wild (WQ) and reintroduced (RQ) individuals. The relative abundance of Actinobacteria was higher in the rhizosphere surrounding wild plants than from reintroduced individuals in Manglietiastrum sinicum [43]. The low abundance of Actinobacteria in the rhizosphere surrounding these reintroduced plants may be a factor contributing to the scarcity of both M. sinicum and P. squamata.

Bradyrhizobium and Streptomyces are known to promote plant growth [44,45]. Bacillus and Pseudomonas control plant disease through having antimicrobial activities [46]. In our study, the relative abundance of Acidothermus was lower in the rhizosphere of the reintroduced P. squamata individuals at Qiaojia (RQ) than in that of the wild individuals at Qiaojia (WQ). We found that the relative abundance of Bacillus was lower in the rhizosphere of wild P. squamata individuals than in that of ex-situ or reintroduced individuals at Qiaojia, which may represent a biotic factor contributing to the scarcity of P. squamata. These PGPR are likely to have crucial effects on the rhizosphere soil bacterial community structure of P. squamata and may have important implications for plant health and survival.

4.3. Relationships between soil physicochemical properties and bacterial community

Soil is crucial in the exchange of organic matter and energy in the soil-microbe-plant ecosystem. Moreover, nutrients can affect disease tolerance or resistances of plants to pathogens [47]. The Bacteroidetes play an important role in the decomposition of polysaccharide organic matter [48]. In our study, we found that both OM content and the relative abundance of the Bacteroidetes were lowest in the EY site, indicating that the soils at the EY site had a low ability to degrade organic matter. Lower OM content tends to indicate low soil fertility, and higher concentrations of disease and pathogens, therefore, in conservation sites, we should regulate and manage soil nutrients to reflect the optimum conditions for P. squamata growth. Because the ex-situ Kunming (EK) had low levels of TK and AN, similar to those in the ex-situ Qiaojia (EQ), this suggests that EK could be a good area for the ex-situ conservation of P. squamata.

4.4. Prediction of function of the P. squamata rhizosphere bacterial community

The functions of the P. dabeshanensis rhizosphere bacteria were mainly related to “amino acid transport and metabolism”, “cell wall/membrane/membrane biogenesis”, “energy production and conversion”, and “signal transduction mechanisms”. “Chemoheterotrophy”, “aerobic_chemoheterotrophy”, “nitrogen fixation”, “cellulolysis” were found to be most abundant [39], which are consistent with our results from P. squamata.

5. Conclusions

This study represents the first exploration of the diversity, composition, and potential function of rhizosphere soil bacterial communities in wild, ex-situ, and reintroduced P. squamata at different conservation sites. Soil TK and AN significantly affected rhizosphere bacterial community structures. However, because the individual trees varied in bacterial community, future studies should include broader sampling of P. squamata individuals for more detailed comparative analysis. A combination of culture-dependent methods to increase levels of plant growth-promoting rhizobacteria (PGPR) and meta-omics to investigate rhizosphere microbiome are potential future conservation tools.