Combined analysis of metagenome and transcriptome revealed the adaptive mechanism of different golden Camellia species in karst regions

Camellia sect. Chrysantha is an important rare and protected plant species. Some golden Camellia species grow in karst soil while others grow in acidic soil. In order to study the adaptation mechanism of golden Camellia to the karst environment, four species of golden Camellia growing in the karst soil (Camellia pubipetala, Camellia perpetua, Camellia grandis, and Camellia limonia) and four species growing in the acidic soil (Camellia nitidissima, Camellia euphlebia, Camellia tunghinensis, and Camellia parvipetala) were selected for this study. Combining the metagenome and transcriptome, the structure and function of the rhizosphere microbial communities and the gene expression in roots of golden Camellia were analyzed. The results showed that the rhizosphere microbial communities in different golden Camellia were significantly different in abundance of Acidobacteria, Actinobacteria, Candidatus_Rokubacteria, Nitrospirae, Planctomycetes, and Candidatus_Tectomicrobia. The proportion of Candidatus_Rokubacteria was significantly higher in the rhizosphere soil of four species of golden Camellia grown in karst areas, compared to C. nitidissima, C. euphlebia, and C. tunghinensis. The linear discriminant analysis Effect Size showed that C. parvipetala was similar to karst species in the enrichment of ABC transporters and quorum sensing. During the transcriptome analysis, numerous upregulated genes in four karst species, including CYP81E, CHS, F3H, C12RT1, NAS, and CAD, were found to be enriched in the secondary metabolite synthesis pathway in the KEGG library, when compared to C. tunghinensis. This study provides information for plant adaptation mechanisms on the rhizosphere soil microbial composition and gene expression in secondary metabolic pathways to karst habitats and its distribution in karst areas.


Introduction
Camellia sect.Chrysantha Chang is an evergreen shrub or small tree of the genus Camellia and the family Theaceae (Wei et al., 2017).They are second-class protected wild plants in China.At present, there are more than 20 species of plants in the Camellia sect.Chrysantha reported, and it is famous for its golden Camellia flowers.According to the different ecological environments and growth soils of golden Camellia, it is mainly divided into acid soil and calcareous soil golden Camellia (Su and Mo, 1988).Acid soil golden Camellia mainly grows in acidic soil, the soil property is mainly acid red soil, and calcareous soil golden Camellia mainly grows in karst soil.Different from the acidic soil, the soil in karst regions is alkaline, with high content of calcium and magnesium, low soil content, and less water (Zhu et al., 2022a).Therefore, high Ca 2+ and drought stress become the main factors limiting plant growth in karst regions.So far, no species of golden Camellia has been found to be able to grow in both acidic and calcareous soil areas.Thus, the adaptation mechanism of golden Camellia to the environment may be species-specific, leading to different environmental adaptation mechanisms among different species of golden Camellia.Therefore, it is necessary to compare the species growing in acidic soil with the ones growing in karst soil to investigate the adaptation mechanism of golden Camellia.
At present, some studies have investigated the adaptation of plants in karst regions, and offered some basic information for karst ecosystem and vegetation protection.Firstly, the water stress in karst regions is an important factor for plant growth.It has been found that the physiological morphology of plants has an important drought resistance.Plants in karst regions often have a thick cuticle by closing epidermal stomata, reducing water evaporation and transpiration (Deng et al., 2004;Rong et al., 2005).Physiologically, karst plants can alleviate the damage caused by drought stress through the accumulation of osmoregulatory substances, including soluble protein, proline, and flavonoids, to enhance their drought resistance (Guo et al., 2018).Secondly, in addition to drought, high calcium is also one of the significant characteristics of soil in karst regions.As an essential element for plant growth and development, calcium plays an important role in signal transduction, but high content of calcium in soil will cause osmotic damage to plant growth (Wang et al., 2011;Li et al., 2014).In order to adapt to the high calcium environment in karst regions, plants have strong regulation of low intracellular calcium content.It was found that plants in karst regions could reduce the effects of Ca 2 + on their growth and development by transporting Ca 2+ to other parts of the body (Jarrett et al., 1982;Dolman et al., 2003;Zhang et al., 2017;Jin et al., 2018).
Certain types of golden Camellia growing in karst regions exhibit some adaptive features to their karst environment.In our previous study, golden Camellia from different growth areas (Camellia nitidissima from acidic soil and Camellia limonia from karst soil) were treated with high calcium to analyze the adaptation mechanism of golden Camellia in karst environment according to their differences in morphology and physiology, and it was found that golden Camellia growing in the karst soil has a more stable photosynthetic capacity and water retention ability (Liu et al., 2022a).Various studies have recently demonstrated that golden Camellia was capable of adapting to drought stress in karst areas by changing the size of its stomata and by developing leathery leaves (Chai et al., 2021;Zhu et al., 2021;Zhu et al., 2022b), but the specific adaptation mechanism in transcription level is still unclear.Moreover, soil microorganisms in different regions are also crucial to plant growth, and there are few studies that focus on the roles of rhizosphere microorganisms in the adaptation mechanism of karst plants.
In this study, a total of eight species of golden Camellia (four in acidic soil and four in karst soil) were used.The structural and functional characteristics of rhizosphere microorganisms in eight different species of golden Camellia growth regions were analyzed by metagenomic analysis, and the roles of rhizosphere microorganisms in the adaptation mechanism of golden camellia were investigated.In addition, regulator genes in plant roots were analyzed by transcriptome.The combination of metagenome and transcriptome was performed to investigate the adaptive mechanisms of golden Camellia to karst regions.This result provides valuable insights into the molecular level of adaptive mechanisms and gene functions in karst plants.

Plant and rhizosphere soil sample collection
In this study, the root and rhizosphere soil samples of eight species of golden Camellia that grow in the karst soil or acidic soil were used as experimental materials in different nature reserves in Guangxi Zhuang Autonomous Region, China.Four species growing in acidic soil were C. nitidissima (Cni), Camellia euphlebia (Ceu), Camellia tunghinensis (Ctu), and Camellia parvipetala (Cpa), and four species growing in karst soil were Camellia pubipetala (Cpu), Camellia perpetua (Cpe), Camellia grandis (Cgr), and Camellia limonia (Cli) (details in Table 1; Figure 1).Three biological replicates were taken for the rhizosphere soil and root samples of each golden Camellia species.A total of 24 rhizosphere soil and 24 root samples were collected (Figure S1).The rhizosphere soil collection included washing the roots with distilled water three to six times to completely remove non-rhizosphere soil.About 0.5 g of rhizosphere soil, which refers to the soil within a few millimeters of the root surface of the golden Camellia, was collected and packed in sterile centrifuge tubes.Root samples were taken from the root tip, meristematic, elongation, and maturation zones, respectively.Prior to DNA and RNA extraction, the samples were frozen in liquid nitrogen and transported to the laboratory on dry ice.The roots of golden Camellia were ultrasonically washed twice with sterile water for 5 min each time to remove soil and microorganisms from the root surface.After washing, approximately 0.5 g of plant root tissue was selected as root samples and collected in a sterile 1.5-mL centrifuge tube for RNA extraction.The soil pH and calcium concentration were determined and published as part of the results (Supplementary Table S1) (Zhu et al., 2022a).The E.Z.N.A. ® Soil DNA Kit from Omega Bio-Tek (GA, USA) was utilized to extract DNA from soil samples.Following DNA extraction, the purity and concentration of the DNA and the length of extracted DNA were measured with NanoDrop2000 from Thermo Fisher Scientific (MA, USA) and using 1% agarose gel electrophoresis.DNA fragmentation was accomplished using Covaris M220 from Shanghai Pudi Biotechnology Co., Ltd.(Shanghai, China).Primers 515F (5′-barcode-GTGCCAGCMGCCGCGCGGGTAA-3′)/806R (5′-GGACTACHVGGGTWTCTAAT-3′) and ITS1F (5′-barcode- ) were used to amplify the V3-V4 region of the bacterial 16S rRNA gene and the ITS1 gene region of fungi, respectively.Amplicons were detected on 2% (w/v) agarose gels and purified using the AxyPrep DNA Gel Extraction Kit (Axygen Biosciences).PE libraries were constructed using the NEXTFLEX ® Rapid DNA-Seq Kit (Bioo Scientific, Texas, USA).The construction process included the following: (1) connecting the joints, (2) removing the joint self-connecting segments using magnetic bead screening, (3) enriching the library templates by PCR amplification, and (4) recovering the library through magnetic beads from the PCR product.The genome sequencing was performed with Illumina Novaseq 6000 (Illumina, CA, USA) in Shanghai Megi Biomedical Technology Co.

Sequence pre-processing, assembly, gene set construction, and abundance calculation
Fastp version 0.23.0 was used for data quality control based on the original sequencing data.Low-quality and N-containing reads were removed to obtain high-quality sequences needed for subsequent analysis.After quality control, the short segment sequence was assembled, and the open reading frames were predicted.After quality control, Megahit was used to assemble multiple short sequences that were mixed together.(Megahit can be found at https://github.com/voutcn/megahit).Sequences with gene sequence length >100 bp were predicted for ORFs using Prodigal (http://compbio.ornl.gov/prodigal).The predicted gene sequences of the samples were clustered, the nonredundant gene set was constructed, and the base sequence of the nonredundant gene set genes was obtained using CD-HIT (http:// www.bioinformatics.org/cd-hit/).SOAPaligner (http:// soap.genomics.org.cn/) was used to assess the quality of reads for each sample, compare the redundant gene sets, and determine the abundance of genes present in the statistical sample.Gene abundance was quantified as Reads Per Kilobase Million (RPKM).

Annotations on species and functions
Non-redundant gene sets from various databases, including NR database, carbohydrate-active enzyme database, and Kyoto Encyclopedia of Genes and Genomes (http://www.genome.jp/kegg/) database, were compared using DIAMOND (https:// github.com/bbuchfink/diamond)with the blastp parameter and an E-value of ≤1e-5.

Transcriptome sequencing and sequence analysis 2.3.1 RNA extraction, library construction and sequencing
Total RNA was extracted from the root samples using Plant RNA Purification Reagent (Invitrogen, Carlsbad, CA, USA).RNA concentration and purity were detected using Nanodrop2000 (Thermo Fisher Scientific, MA, USA), while RNA integrity was verified via 1% agarose gel electrophoresis.Then, 5 mg of total RNA from each sample was used for library preparation with a TruSeq Stranded Total RNA Sample Preparation Kit (Illumina, San Diego, USA).RIN value was determined using Agilent 2100 (Agilent Technologies, Palo Alto, CA, USA).The fragment RNA was then reversely transcribed into cDNA, then connected to adaptor, and the original data were sequenced on Illumina Novaseq 6000 (Illumina, CA, USA).

Data preprocessing, comparison, and expression analysis
Clean reads were obtained by removing raw reads from the original data for subsequent analysis.Clean data (reads) after quality control were used for TopHat2 (http://tophat.cbcb.umd.edu/)software comparison with reference genome (Camellia sinensis, https://www.ncbi.nih.gov/genome/?term=GCF_004153795.1) to get mapped data (reads).Cufflinks (http://cole-trapnelllab.Making.IO/ Cufflinks/) was used to assemble the mapped reads splicing.After genes were identified in the gene ontology (GO) and KEGG database, the information of the function of genes and transcription was obtained.RSEM (http://deweylab.biostat.wisc.edu/rsem/)was used to quantitatively analyze the gene expression levels, and the quantitative index was the transcripts per million reads (TPM).DESeq2 was used for inter-sample or inter-group gene differential expression analysis, and the differentially expressed genes (DEGs) were identified.In this study, genes and transcripts with absolute multiple change ≥2 and false discovery rate (FDR) < 0.05 were DEGs.

Statistical analysis
TBtools (version 1.047) was used for heat map analysis of DEGs.The Wilcoxon rank sum test was used to detect significant differences in gene abundance in soil samples in metagenome data, and the mean was analyzed at the p < 0.05 significance level.

Metagenome and transcriptome datasets of eight species of golden Camellia growing in karst soil and acidic soil
After quality control, the metagenomic data yielded approximately 116 million clean reads, averaging approximately 5 million clean reads per sample.In each sample, clean reads comprised over 98% of the raw reads.About 1 million contigs were obtained with an average of 443,771 contigs per sample.The high quality of sequence assembly was reflected in N50 and N90 values.Prediction of open reading frames (ORFs) yielded about 1.2 million sequences, ranging from 102 to 9,627 bp.The characteristics of the DNA-sequencing results indicate that the data are of high quality and suitable for further analysis (Supplementary Table S2).
Similarly, for the transcriptome data after quality control, a total of 165.05 Gb clean data was obtained, and the clean data of each sample reached more than 5.87 Gb, the percentage of Q30 base was more than 91.64%, and the error rate of each sample was about 0.2%.GC content ranged from 44.08% to 47.64%, which met the requirements of subsequent analysis.The clean reads of each sample were aligned with the designated reference genome, and the alignment rate ranged from 68.02% to 76.58% (Supplementary Table S3).At the same time, the gene expression level of each sample was in line with the sequencing precision range of transcriptome, and the data were evenly distributed (Figure S2).

Characteristics of rhizosphere microbial communities of golden Camellia growing in the acidic soil and karst soil
There are many important factors that affect the growth of golden Camellia in specific habitats (Figure 1).Not only in terms of phenotype, but also the characteristics of soil microbial community are closely related to the habitat adaptation of golden Camellia.Based on metagenome analysis, the characteristics of rhizosphere microbial community structure of different golden Camellia growing in the acidic soil or karst soil were examined, and their roles were analyzed.Figure 2A shows that the distances of the three biological replicates of the same golden Camellia were close, indicating the reliability of metagenome sequencing data.The soil microbial community structure of golden Camellia (Cni, Ceu, and Ctu) in the same habitat (acidic soil) has a similarity.Moreover, the rhizosphere microbial community structure of golden Camellia (Cpu, Cpe, Cgr, and Cli) growing in karst soil was also similar.It is worth noting that Cpa growing in the acidic soil was closer to the karst soil in terms of community structure.
Based on the genes, the species classification annotation of eight golden Camellia species of rhizosphere microorganisms was statistically analyzed (Table 2).In acidic soil (Cni, Ceu, Ctu, and Cpa), the number of microbial phyla ranges from 200 to 209, and the number of microbial species ranges from 16,714 to 18,208.In karst soil (Cpu, Cpe, Cgr, and Cli), the number of microbial phyla varied between 208 and 217 and the number of microbial species varied between 19,454 and 21,917.These results demonstrated that the species diversity of rhizosphere microorganisms in different habitats was different, and the diversity of rhizosphere microorganisms in karst soil was higher than that in acidic soil.The rhizosphere microbial composition of eight golden Camellia species based on phylum and class levels was analyzed (Figures 2B-D).At the phylum level, the top three phyla with higher abundance for eight golden Camellia species were all Proteobacteria, Actinobacteria, and Acidobacteria.However, the mean proportions of microbial phyla in different golden Camellia were highly significantly different, particularly in Acidobacteria, Actinobacteria, Candidatus_Rokubacteria, Nitrospirae, Planctomycetes, and Candidatus_Tectomicrobia (0.001< p-value ≤ 0.01) (Figure 2E).Among the top five microbial phyla in abundance, only Acidobacteria was significantly more abundant in normal species (Cni, Ceu, and Ctu) than in Cpa, Cpu, Cpe, Cgr, and Cli.The abundance of Candidatus_Rokubacteria was significantly higher in Cpa, Cpu, Cpe, Cgr, and Cli than in Cni, Ceu, and Ctu.This result indicated that mainly Acidobacteria and Candidatus_Rokubacteria have significant abundance differences between rhizosphere soil of normal and karst soils.Furthermore, the mean proportions of Acidobacteria and Actinobacteria also significantly differed in Cni, Ceu, and Ctu growing in acidic soil, among which the mean proportion of Acidobacteria in Cni was the highest (47.36%), and the mean proportion of Acidobacteria in Ceu was the lowest (23.56%).Among Actinobacteria, Ceu was the highest (17.41%) and Cni was the lowest (8.214%) (Figure S3A).Nevertheless, there were no significant differences in microbial phyla in karst soil, except for Nitrospirae which had low mean proportions (Figure S3C).
At the class level, Alphaproteobacteria was the most dominant bacteria in the eight golden Camellia species (Figures 2B-D).Besides Alphaproteobacteria, the dominant bacteria include unclassified_p_Acidobacteria, Actinomycetia, Acidobacteriia and Gammaproteobacterial in Cni, Ceu, and Ctu.Unclassified_p_ Candidatus_Rokubacteria, unclassified_p_Acidobacteria, Deltaproteobacteria, Betaproteobacteria, and Actinomycetia were dominant bacterial community in Cpa.Unclassified_p_ Candidatus_Rokubacteria, Actinomycetia, Deltaproteobacteria, unclassified_p_Acidobacteria, unclassified_p_Actinobacteria, Betaproteobacteria, and Thermoleophilia were dominant bacterial community in karst species (Cpu, Cpe, Cgr, and Cli).The mean proportions of Alphaproteobacteria, Deltaproteobacteria, Betaproteobacteria, and Gammaproteobacterial differed significantly in eight golden Camellia species (Figure 2F).The mean proportions of Alphaproteobacterial ranged from 22.77% (Cni) to 40.49% (Ceu) in acidic soil, and it ranged from 12.3% (Cpu) to 24.72% (Cli) in karst soil.However, the proportions of Deltaproteobacteria and Betaproteobacteria in karst soil were significantly higher than those in acidic soil (Cni, Ceu, Ctu).Only the proportion of Alphaproteobacteria was different in the same habitat (karst soil) (Figures S3B, D).At phylum and class levels, the differences in microbial community structure between the same or different habitats may reflect the functional characteristics of soil microorganisms.

Functional pathways of rhizosphere microorganisms of golden Camellia growing in the normal and karst soils
The main functional pathways of rhizosphere microorganisms could have a certain impact on plant growth.Therefore, the functional pathways of soil microorganisms were analyzed.The functional annotation of non-redundant genes in all samples was performed to compare the abundance of functional genes in the rhizosphere microorganisms of eight golden Camellia species based on the Kyoto Encyclopedia of Genes and Genomes (KEGG) database and carbohydrate-active enzyme (CAZy) functional database.In total, there were 111,555,026 KEGG pathway-associated genes (Supplementary Table S4) and 6,721,670 CAZy genes (Supplementary Table S5) from all rhizosphere soil samples.The global and overview maps, the carbohydrate metabolism, amino acid metabolism, and energy metabolism were the main functional pathways based on pathway level 2 (Figure 3A).Based on the KEGG pathway level 3, all soil samples were divided into three main types (1: Cni and Ctu; 2: Ceu, Cpe, and Cli; 3: Cpa, Cpu, and Cgr) (Figure 3B).Those were divided into two types based on CAZy (1: Cni, Ceu, and Ctu; 2: Cpa, Cpu, Cpe, Cgr, and Cli) (Figure 3C).This reflected that the functions of soil microorganisms of the same type were somewhat similar and Cpa was special with more similarity with karst species.
The LEfSe analysis was performed on all KEGG pathway level 3 to identify pathways with significant differences in non-redundant gene sets after KEGG annotation (Figures 3D-H).In normal species (Cni, Ceu, Ctu, and Cpa), Cni was mainly related to the metabolic pathway, while Cpa was mainly related to ABC transporters, quorum sensing, and biosynthesis of amino acids.However, no significantly different pathways were identified in Ceu and Ctu (Figure 3D).In addition, in karst species (Cpu, Cpe, Cgr, and Cli), Cpu was mainly related to biosynthesis of secondary metabolites.
Cli was mainly related to quorum sensing, microbial metabolism in diverse environments, two-component system, and ABC transporters.Similarly, no significantly different pathways were identified between Cpe and Cgr (Figure 3E).When Cni (normal species) was compared with Cpu, Cpe, Cgr, and Cli, respectively, Cni was associated with the metabolic pathway in each comparison.

The correlations between rhizosphere microbial community structure and function in normal and karst soils
The function of rhizosphere microorganisms under different golden Camellia is related to its composition.Species and functional contribution analysis further confirmed the relationship between rhizosphere microorganisms at the class level and KEGG pathway level 3 (Figure 3J).In normal species Cni, Ceu, and Ctu, the top four contributing microbial classes in the metabolic pathway were Alphaproteobacterial (43.56%-28.47%),Actinomycetia (17.34%-9.11%),unclassified_p_Acidobacteria (20.26%-9.23%),and Acidobacteriia (17.93%-9.04%).Therefore, the high abundance of Alphaproteobacteria, unclassified_p_Acidobacteria, and Acidobacteriia in Cni, Ceu, and Ctu resulted in a higher abundance of metabolic pathway than Cpa, Cpu, Cpe, Cgr, and Cli.Quorum sensing and ABC transporters were mainly related to rhizosphere microorganisms in karst soil.Alphaproteobacteria (33.84%-22.38%),Actinomycetia (12.70%-7.03%),and unclassified_p_Candidatus_Rokubacteria (31.02%-23.77%)were the major contributing classes to quorum sensing.The ABC t r a n s p o r t e r s p a t h w a y w a s m a i n l y a s s o c i a t e d w i t h Alphaproteobacteria (37.19%-25.37%),unclassified_p_ C a n d i d a t u s _ R o k u b a c t e r i a ( 2 2 .3 3 % -1 5 .9 6 % ) , a n d Deltaproteobacteria (14.93%-9.84%).Overall, the functional contribution patterns of soil microorganisms in Cpa were similar to those in golden Camellia growing in karst soil.

The correlations of root gene expression among the eight species of golden Camellia
The correlation of root gene expression of eight species of golden Camellia was performed (Figure 4A).Cpa was semblable to that of Cpu (0.749) and Cgr (0.759).There was comparability in Cni, Ceu, and Ctu because the correlations among Cni, Ceu, and Ctu were all above 0.7.Furthermore, the correlation among golden Camellia growing in karst soil was also above 0.678.In particular, Cli and Cgr showed a high correlation of 0.889.
In Cni, Ceu, Ctu, and Cpa, there were 18,236 (54.99%) common genes, and Cpa had the highest number of unique genes (2,372),  Liu et al. 10.3389/fpls.2023.1180472Frontiers in Plant Science frontiersin.orgwhich was consistent with the metagenomic results (Figure 4B).Among Cpu, Cpe, Cgr, and Cli growing in karst soil, there were 20,420 (62.30%) common genes, which might be related to the adaptation of golden Camellia to karst habitats (Figure 4C).

DEGs of root gene expression among different golden Camellia
Cni, Ceu, and Ctu were selected and compared with golden Camellia (Cpu, Cpe, Cgr, Cli, and Cpa), respectively, to investigate the differences in the adaptation of different golden Camellia to karst habitats (Figure 5).In the comparisons of DEGs in Cpa_vs_Cni, Cpu_vs_Cni, Cpe_vs_Cni, Cgr_vs_Cni, Cli_vs_Cni, Cpa_vs_Ceu, Cpu_vs_Ceu, and Cpe_vs_Ceu, the numbers of downregulated genes were higher than the numbers of upregulated genes in each comparison group (Figures 5A, B; Table S1).

The annotation and functional enrichment of DEGs in different groups
The annotation and functional enrichment of DEGs were performed by the GO and KEGG database to further analyze the regulation pathways of golden Camellia in different habitats.Based on the above DEGs, the largest number of DEGs in the Ctu-related group was selected for annotation and functional enrichment in this study (Figures 6, S5).
In GO and KEGG database annotation, the distribution trend of DEGs in each Ctu-related comparison group was basically the same in each GO term and KEGG functional pathway (Figures 6A, B 6C-G).In Cpa_vs_Ctu, non-homologous end-joining was the most enriched (rich factor = 0.43), while nucleocytoplasmic transport was the most enriched (number = 106) (Figure 6C).Also, non-homologous endjoining was also the most enriched in Cpu_vs_Ctu (rich factor = 0.52) (Figure 6D).Among DEGs in Cpu_vs_Ctu, they were enriched in phenylpropanoid biosynthesis, stilbenoid, diarylheptanoid and gingerol biosynthesis, and flavonoid biosynthesis way related to secondary metabolites, etc. Pathways related to secondary metabolites were enriched not only in DEGs of Cpu_vs_Ctu, but also in those of Cpe_vs_Ctu, Cgr_vs_Ctu, and Cli_vs_Ctu, especially in phenylpropanoid biosynthesis (number: Cpu_vs_Ctu = 148; Cpe_vs_Ctu = 142; Cgr_vs_Ctu = 128; Cli_vs_Ctu = 99) (Figures 6D-G).These results suggested that the adaptation of golden Camellia growing in karst soil is one of the mechanisms by regulating its own metabolic substances.

The secondary metabolic biosynthetic pathways in different golden Camellia species growing in karst soil
After annotation and functional enrichment of the DEGs in the KEGG database, the results showed that the DEGs were significantly enriched in secondary metabolic biosynthetic pathways.Gene expression profiles were shown in the pathway of secondary metabolite synthesis to illustrate the adaptation of golden Camellia to karst habitats, which was done using Ctu as the control (Figure 7).The DEGs in the Ctu-related comparison groups were mainly upregulated in flavonoid biosynthesis, phenylpropanoid biosynthesis, isoflavonoid biosynthesis, and nitrogen metabolism.In the isoflavonoid biosynthesis, there were 11 DEGs related to CYP81E enriched, namely, 8 upregulated genes and 3 downregulated genes.In flavonoid biosynthesis, the upregulated genes include CHS and CHI.In addition, the upregulation of genes involved in flavonoid biosynthesis accelerated the production of various flavonoid species, such as C12RT1, which was upregulated on Cpa, Cpu, Cpe, Cgr, and Cli.Compared to Ctu, the expressions of CAD and POD in Cpa, Cpu, Cpe, and Cli were increased.In soil rhizosphere microorganisms, there were also related DEGs in the process of nitrogen cycle, and the upregulated genes were more than the downregulated genes, such as amoA and nirB.This suggests that rhizosphere microorganism may assist the adaptation of golden Camellia to karst habitats by accelerating nitrogen cycling.

Discussion
By combining metagenome and transcriptome approaches, we analyzed the differences in rhizosphere microbial structure and gene expression in roots of golden Camellias in order to understand the adaptation mechanisms of golden Camellia in karst environments.
First of all, in terms of rhizosphere microbial composition, Cni, Ceu, and Ctu growing in acidic soil have high similarity, while Cpa has more similarity with Cpu, Cpe, Cgr, and Cli, which grow in karst soil.In addition, Cpa is the more drought-resistant golden Camellia growing in acidic soil (Su and Mo, 1988).In terms of rhizosphere soil microbial composition, the proportion of Proteobacteria at the phylum level was the highest among the eight species of golden Camellia.Proteobacteria is the largest bacterial phylum, including Alphaproteobacteria, Betaproteobacteria, Deltaproteobacteria, and Gammaproteobacteria, which can survive under stress conditions such as drought, heavy metals, and acidic environment (Wang et al., 2019;Liu et al., 2021).The results of this study also found that Proteobacteria, especially Alphaproteobacteria, had the highest functional contribution in the top seven abundance in rhizosphere soil (Figure 3J).Alphaproteobacteria include nitrogen-fixing and nitrifying bacteria that reduce ammonia and ammonium to nitrate, while Betaproteobacteria reduce nitrate to nitrite (Einsle, 2011;Li et al., 2012).Ceu had the highest proportion of Alphaproteobacteria, followed by Ctu, and genes related to NH 4 + reduction to NO 3 − , such as hao and nxrA, also had higher gene abundance in the Cturelated comparison groups (Figures 2F, 7), while the proportion of Betaproteobacteria was higher in Cpa and Cpu.The abundance of Alphaproteobacteria had been reported to be inversely proportional to high nitrogen concentration (Li et al.,FIGURE 7 Gene expression of golden Camellia growing in karst soil in the secondary metabolite synthesis pathways. 2019; Zhu et al., 2022a).Moreover, at the class level, Deltaproteobacteria and Betaproteobacteria were more abundant in karst soil than in acidic soil.Acidobacteria belonged to acidophilus, and its proportion in Cni, Ceu, and Ctu was significantly higher than that in golden Camellia growing in karst soil.As an important component of microbial community, it has important ecological functions by promoting carbon, nitrogen, and sulfur metabolism in soil (Kalam et al., 2020).Different from Acidobacteria, the proportion of Candidatus_Rokubacteria in Cpa, Cpu, Cpe, Cgr, and Cli was significantly higher than that of Cni, Ceu, and Ctu.Candidatus_Rokubacteria growth and reproduction are promoted by soil nitrogen, calcium, magnesium, and other elements (Ivanova et al., 2021;Zhu et al., 2022a).Among the functions of rhizosphere soil microorganisms in karst habitats, ABC transporters and quorum sensing have scores in high linear discriminant analysis, and the two most important contributors to these functions were Alphaproteobacteria and Candidatus_Rokubacteria (Figure 3).These results indicated that the karst habitats affected the survival of Candidatus_Rokubacteria and Candidatus_Rokubacteria significantly affected the function of soil microbial community in karst habitats.Actinobacteria was the third most abundant phylum at the phylum level, with higher proportions in Cpu, Cpe, and Cgr (Figure 2E).Actinobacteria have the ability to decompose organic matter, especially cellulose and chitin (Liu et al., 2014;Liu et al., 2016;Li et al., 2019;Zhu et al., 2022a).In addition, the decomposition of nitride-containing substances can increase the nitrogen content in soil, which may enhance the absorption of carbon and nitrogen by plants, and the absorption of nitrogen in plants can promote the production of secondary metabolites, thereby alleviating the damage caused by external stress (Figure 8).
Soil properties affect the composition of rhizosphere soil microorganisms in different habitats, which, in turn, affect their functional pathways and may assist in the adaptation of golden Camellia to karst habitats.The results showed that ABC transporters and quorum sensing were most different in rhizosphere soil microorganisms of golden Camellia growing in karst soil compared with those growing in acidic soil.As one of the largest protein families, ABC transporters are responsible for transporting carbohydrates, lipids, and proteins (Qu et al., 2020).Significantly, ABC transporters in rhizosphere microorganisms can mediate the exchange of substances between the roots and soil, which may be beneficial to the absorption of ferrite, phosphorus, and zinc ion by the roots of plants, thus promoting its growth and development (Ryan et al., 2001).Quorum sensing is due to a certain density of microorganisms, which allows bacteria to communicate with each other (Basavaraju et al., 2016).CAZy genes include glycosyl transferases (GT), glycoside hydrolases (GH), carbohydrate esterases (CE), auxiliary activities (AA), polysaccharide lyases (PL), and carbohydrate-binding modules (CBM), which have the functions of degradation and modification.GT and GH accounted for the highest proportion in the rhizosphere soil microorganisms of the eight species of golden Camellia.GT is responsible for the synthesis of sugars and sugar-containing substances (Lairson et al., 2008).GH has the ability to decompose organic matter to form various monosaccharides and oligosaccharides, which provide energy for microorganisms themselves and plants (Manoharan et al., 2015).
When plants are exposed to abiotic stress, their gene expression levels will also change in response to external stress.Drought and high calcium stress in karst regions have a certain negative effect on the growth of plants, which are not adaptive to karst regions.Calcium is an essential mineral element for plant growth and development.It plays an important role in cell components, signal transduction, and other aspects (Yang et al., 2008).In addition, when plants are exposed to high calcium environment for a long time, the synthesis of various including organic acids and osmotic substances, helps plants adapt to the high calcium environment (Zhang et al., 2021).Studies have shown that the secondary metabolic contents of flavonoids, anthocyanins, saponins in apple, and Panax notoginseng are increased in a high-calcium environment (Zu et al., 2022).High concentration of calcium in the soil will cause an osmotic toxic effect on the growth of golden Camellia, reducing the activities of various antioxidant enzymes, such as superoxide dismutase (SOD), catalase (CAT), and peroxidase (POD), causing serious oxidative damage to plants (Yang and Fang, 2022).In the present study, KEGG enrichment analysis of DEGs in the Ctu-related comparison groups, especially in Cpu_vs_Ctu, Cpe_vs_Ctu, Cgr_vs_Ctu, and Cli_vs_Ctu, was found to be enriched in multiple secondary metabolic synthesis pathways, including phenylpropanoid biosynthesis, flavonoid biosynthesis, and isoflavonoid biosynthesis.Phenylpropanoid, flavonoid, and isoflavonoid can eliminate free radicals in plants and enhance their antioxidant capacity to reduce oxidative damage caused by high calcium environment (Li et al., 2018).Moreover, some flavonoids, such as anthocyanins, can enhance the color of plant leaves and maintain a high photosynthetic capacity of plants, which is an important adaptive mechanism for golden Camellia to adapt to the high-calcium environment (Liu et al., 2022a).Firstly, in isoflavonoid biosynthesis, eight genes were found to be upregulated and related to CYP81E, which is an important metabolic enzyme in various isoflavones' synthesis in roots of Cpu, Cpe, Cgr, and Cli compared with Ctu.The expression of CYP81E-related genes was mainly related to the accumulation of glyceollin I/II/III, 2′-Hydroxygenistein, 2′-Hydroxybiochanin A, vesritol, medicarpin-3-O-glucoside-6-Omalonate, and pisatin (Liu et al., 2003;Liu et al., 2022b;Wang et al., 2023).Likely, during soybean germination, glyceollin I accumulation can eliminate reactive oxygen species (ROS) caused by microbial stress and alleviate membrane lipid peroxidation (Jeon et al., 2012).The synthesis of pisatin can enhance the resistance of pea seeds to osmotic stress (Brosowska-Arendt et al., 2014).The expression of CHS in flavonoid biosynthesis pathway was downregulated in Feizixiao lychee fruit after high calcium treatment which may cause oxidative damage (Shui et al., 2022).Therefore, the upregulated expression of CHS in the roots of golden Camellia growing in karst soil may help its adaptation to the high Ca 2 + environment.The synthesis of a variety of flavonoids may promote the adaptation of golden Camellia to karst habitats.Flavonoids and isoflavones have similar effects on plants and can mitigate oxidative damage caused by scavenging ROS produced when plants are stressed (Hernandez et al., 2009).In drought-tolerant rice cultivars, oxidative damage induced by drought stress is alleviated by the accumulation of flavonoids, including naringin (Dey and Bhattacharjee, 2020).In the roots of salt-tolerant wheat varieties, genes related to phenylpropanoid biosynthesis, including 4CL and CAD, were significantly upregulated (Yue et al., 2022).In the phenylpropanoid biosynthesis pathway, the upregulated expression of CAD and POD involved the production of various lignin.Lignin can provide enhanced cell support that improves plant body support and water removal, thereby increasing plant drought resistance (Cabane et al., 2012).A complex interplay exists between the nitrogen metabolic pathway and secondary metabolite synthesis.The synthesis of secondary metabolites such as flavonoids, isoflavones, and phenylpropanoids derived from phenylalanine and tyrosine is also affected by the uptake of carbon and nitrogen, so there is a complex interaction between the nitrogen metabolic pathway and the synthesis of secondary metabolites (Garcia-Calderon et al., 2020).The genes related to NRT and NR can enhance the uptake of extracellular nitrate and nitrite, with nitrogen cycling in the rhizosphere soil facilitated by microorganisms (Zhao et al., 2018).In the present study, the roots of golden Camellia growing in karst areas demonstrated the upregulation of genes related to NRT and NR.Golden Camellia's adaptation to karst regions may associate with the rhizosphere microbe composition and the expression of genes in root, especially the genes that were involved in secondary metabolite synthesis pathways.

Conclusion
In this study, metagenomic and transcriptome sequencing were used to analyze the differences in rhizosphere soil microbial composition and root gene expression between golden Camellia growing in acidic soil and that growing in karst soil to compare and investigate their differences in adaptation to karst areas.Compared to the expression level of root genes in golden Camellia growing in acidic soil, those in karst soil are upregulated for genes in the secondary metabolic synthesis pathway and for genes involved in nitrogen transport.Thus, through the absorption of some nitrogen substances to adapt to the karst habitat, golden Camellia may enhance secondary metabolite synthesis.The abundance of Candidatus_Rokubacteria in rhizosphere soil microorganisms of golden Camellia growing in karst soil is significantly higher than that in acidic soil, which makes them have stronger ABC transporters and quorum sensing functions that may affect the transport of materials.Soil microorganisms in the rhizosphere of golden Camellia grown in karst soil exhibit a significant capability to improve soil nitrogen cycling.Rhizosphere microorganisms have the potential to enhance nitrogen absorption in the roots of golden Camellia by improving material transport in the soil.Consequently, this assists golden Camellia in synthesizing secondary metabolites, which could significantly contribute to the survival and growth of golden Camellia in karst regions.manuscript; HJ, YL, and SC also collected the samples.All contributed to the article and approved the submitted version.
FIGURE 2 Characteristics of rhizosphere microbial communities of golden Camellia growing in normal or karst soil.(A) Typing analysis for all soil samples based on the genus level of rhizosphere microorganisms and the distance algorithm is JSD.(B-D) The relative abundance of microorganisms in eight golden Camellia species for major taxa at phylum (upper) and class level (lower).(E, F) The microbial abundance at the phylum and class level in eight golden Camellia species.Star means p-value ≤ 0.05.
FIGURE 3 Functional pathways of rhizosphere microorganisms of golden Camellia growing in the normal and karst soils.(A) KEGG functional annotation for all soil samples was based on pathway levels 1 and 2. The abundance was calculated as Reads Number.(B, C) Typing analysis for rhizosphere microbial gene functions was based on KEGG (B) and CAZy (C) databases.(D-H) LEfSe analysis at KEGG pathway level 3 with significant differences in rhizosphere microbial gene function.Cni, Ceu, Ctu, and Cpa (D).Cpu, Cpe, Cgr, and Cli (E).Cni and Cli (F).Cni and Cpa (G).Cpa and Cli (H).The threshold is 3. (I) The functional composition of rhizosphere microbial community was based on the CAZy database.(J) Functional contribution analysis of the top seven functions at KEGG pathway level 3 based on the abundance of rhizosphere microbial community at the class level.
FIGURE 4 Transcriptome sequencing analysis of root samples of golden Camellia.(A) Correlation analysis of eight golden Camellia species.(B) The Venn diagram of the root genes at Cni, Ceu, Ctu, and Cpa.(C) The Venn diagram of the root genes at Cpu, Cpe, Cgr, and Cli.

FIGURE 8
FIGURE 8Model of effects of soil properties on rhizosphere microbial communities and their role in nitrogen cycling.The blue arrow indicates the benefit of bacterial growth, and the red arrow represents the inhibition of bacterial growth.

TABLE 1
Sampling information of eight species golden Camellia.

TABLE 2
The information of species classification and annotation in rhizosphere microorganisms of 8 species golden Camellia.