Background {#Sec1}
==========

*Passiflora*, a genus within the Passifloraceae family, comprises approximately 500 species, most of which are distributed throughout Latin America; many species of this genus are widely cultivated for their edible value, medicinal efficacy and ornamental properties \[[@CR1]--[@CR3]\]. In the early twentieth century, *Passiflora edulis* Sims (also known as passion fruit), the newest entry on the fruit market at the time, was introduced in southern China, mainly in Taiwan, Guangdong, Guangxi and Fujian. This exotic and perennial vine is adapted to warm and humid climates and is thus generally intolerant of cold temperatures. However, we have developed a cold-tolerant variety of *P. edulis* ('Pingtang 1') that can survive subzero temperatures and is well adapted to the Karst habitats of Guizhou.

Karst topography is characterized as a typical landscape shaped from the dissolution of soluble rocks, such as limestone (L), dolomite (D), and gypsum \[[@CR4]\]. As one of the most fragile terrestrial ecosystems \[[@CR5], [@CR6]\], Karst landforms are especially developed and widespread in soluble rock areas, such as the Nullarbor Plain in Australia, the Western Highland Rim in the United States, and the Yunnan-Guizhou Plateau in China. In Southwest China, Karst rocky desertification has become a serious threat to local agriculture, forestry, and livestock husbandry \[[@CR7], [@CR8]\], and Karst soils have poor water-holding capacity and nutrient cycling ability. The new *P. edulis* variety 'Pingtang 1' is widely cultivated in the Karst areas of Guizhou. Because of breeding and popularization, this variety has undergone several rounds of testing under extreme low temperature conditions and is well adapted to Karst soils despite their poor water-holding capacity and soil nutrient cycling ability \[[@CR9], [@CR10]\]. However, information regarding the regulatory patterns and signaling pathways of genes in this exotic and perennial vining species related to adaptation to Karst soils remains unknown.

Because of the rapid development of next-generation sequencing (NGS) technology and the completion of whole-genome sequencing projects for many plant species \[[@CR11], [@CR12]\], increased availability of genomic and transcriptomic data has provided insights into root-soil interactions. *P. edulis*, also known as passion fruit and native to South America, is now widely cultivated throughout southern China for its edible fruit, medicinal efficacy and ornamental properties. In this study, cuttings of 'Pingtang 1' were cultivated in L rocky desertification areas and sandy D rocky desertification areas. Nutrient elements of both plants (roots, stems and leaves) and soils in the two plots were systematically analyzed. Furthermore, transcriptome libraries from the roots were constructed and sequenced, and 1314 DEGs (531 upregulated and 801 downregulated) were identified as being significantly differentially expressed between the L and D root cDNA libraries (L_R and D_R, respectively). Most differentially expressed genes (DEGs) were involved in nutrient element uptake and utilization and signal regulation. DEGs in KEGG pathway of plant hormone signal transduction, phenylpropanoid biosynthesis, and biosynthesis of unsaturated fatty acids may reveal the potential responses of plant growth and chilling-resistance in two soil sites. The results provide a global overview of gene expression patterns associated with *P. edulis* adaptability to Karst soils and could contribute to revealing the adaptive mechanism at the molecular level.

Results {#Sec2}
=======

Nutrient elements in soils and plants {#Sec3}
-------------------------------------

Analyses of nutrient elements were based on the average values of three biological replications. Yellow soil is characterized as being acidic; the soils of the L and D plots had pH values of 6.45 and 5.75, respectively, both of which are lower than 7.0, indicating the acidic properties of the two soils.

For soil nutrient elements, the amount of available nitrogen (AN) in two sites were similar, exchangeable calcium (ECa), and exchangeable magnesium (EMg) at the L site were greater than those at the D site. However, the available potassium (AK) content exhibited the opposite trend (Fig. [1](#Fig1){ref-type="fig"}a), and the available phosphorus (AP) content was much higher at the D site than at the L site (Fig. [1](#Fig1){ref-type="fig"}a). The difference in EMg between the two soil samples was significant (*p* \< 0.05), and the differences in AP and, AK were highly significant (*p* \< 0.01).Fig. 1**a**) The nutrient element (AN: alkali hydrolyzable nitrogen, AP: available phosphorus, AK: available potassium, ECa: exchangeable calcium, EMg: exchangeable magnesium) concentrations in the soils at two sites (L and D) (mg/kg); **b**) Percentage of nutrient element (N: nitrogen, P: phosphorus, K: potassium, Ca: calcium, Mg: magnesium) concentrations in the roots at the two sites (L and D); **c**) Percentage of nutrient element (N, P, K, Ca and Mg) concentrations in the stems; **d**) Percentage of nutrient element (N, P, K, Ca and Mg) concentrations in the leaves. One asterisk in the figure indicates a significant difference in elements in the roots or other tissues between sites D and L (*p* \< 0.05), while two asterisks indicate a highly significant difference (*p* \< 0.01)

Generally, the proportions of elements in plant tissues indicated that the total amount of N, K, and Ca were greater than P and Mg (Fig. [1](#Fig1){ref-type="fig"}b, c, d). The amounts of N and K were higher in the D site plants than in the L site plants, while the amounts of Ca, and Mg were higher in the L site plants (Fig. [1](#Fig1){ref-type="fig"}b, c, d). The differences in N, K, Ca, and Mg in the roots growing in the two soil types were highly significant (*p* \< 0.01) (Fig. [1](#Fig1){ref-type="fig"}b). The stems from the two soil types presented highly significant (*p* \< 0.01) Mg concentrations, while the K and Ca concentrations were significantly different (*p* \< 0.05) (Fig. [1](#Fig1){ref-type="fig"}c). In the leaves, N was highly significant (*p* \< 0.01), while K, Ca, and Mg were significant (*p* \< 0.05) (Fig. [1](#Fig1){ref-type="fig"}d).

RNA sequencing and transcriptome assembly {#Sec4}
-----------------------------------------

A total of 253,078,350 raw reads and 244,705,162 clean data were obtained from the four cDNA libraries after sequencing. The mean Q30 and GC bias values were 90.31 and 44.92%, respectively, indicating high-quality sequencing (Table [1](#Tab1){ref-type="table"}). After sequence assembly, 142,677 transcripts ranging from 201 to 16,965 bp in length were obtained. The N50 of all transcripts was 2046 bp, while that of 84,198 unigenes was 2326 bp. The proportion of transcripts whose length ranged from 200 to 500 bp accounted for more than half of all the transcripts, and the different length intervals of the unigenes were evenly distributed (Table [2](#Tab2){ref-type="table"}).Table 1Summary of RNA-seq data from four RNA libraries of two-year-old root of 'Pingtang 1' from two different soil sites in Karst habitsLibrariesRaw ReadsClean ReadsClean BasesQ30 (%)GC (%)L_R163,256,21861,695,4169.25G92.9944.74L_R260,417,09458,150,9188.72G87.5345.18D_R166,971,47465,129,3309.77G92.8045.20D_R262,433,56459,729,4988.96G87.9344.57Summary253,078,350244,705,16236.7GTable 2Distribution of length interval of transcripts and unigenes obtained from the RNA-seq of roots of 'Pingtang 1'Length distributionTranscriptsUnigenesNumberPercentageNumberPercentage200--50080,93756.73%23,42827.82%500-1 k22,54515.80%21,58825.64%1-2 k18,54713.00%18,53422.01%\> 2 k20,64814.47%20,64824.52%Total142,67784,198Total nucleotide135,742,661119,695,604Mean length9511422N5020462326

Gene annotation and functional classification {#Sec5}
---------------------------------------------

There were 52,229, 41,877, 41,644, 41,034, 40,640, 21,234, and 15,645 significantly matched unigenes in the NR, SwissProt, GO, NT, Pfam, KO, and KOG databases, respectively, with successful annotation percentages of 62.03, 49.73, 49.45, 48.73, 48.26, 25.21, and 18.58%, respectively. The NR database yielded the greatest proportion of matched unigenes, while the KOG database yielded the lowest. A total of 9093 unigene hits were universally obtained across the seven databases tested (Table [3](#Tab3){ref-type="table"}). Based on the annotations, the unigenes were further classified into different groups, with 41,633, 15,645, and 21,234 unigenes enriching 50, 25, and 19 terms for the GO, KOG, and KEGG databases, respectively (Fig. [2](#Fig2){ref-type="fig"}). GO analysis revealed enrichment in cellular processes, metabolic processes (biological process \[BP\]), binding, catalytic activity (molecular function \[MF\]), and cell and cell parts (cellular component \[CC\]), thus underlying three additional significant terms (single-organism process (BP), organelle (CC) and macromolecular complex). For the KOG analysis, the top three terms were \[O\] posttranslational modification, protein turnover, and chaperones; \[R\] general function prediction only; and \[J\] translation, ribosomal structure and biogenesis. Considering the potential relationship of gene expression in the different soil types, we focused on the fourth term (\[T\] signal transduction mechanisms) for further analysis. The KEGG analysis provided a feasible reference for linking transcriptomes to biological systems by pathway mapping. We focused on the following top five terms: carbohydrate metabolism; translation; overview of metabolism; folding, sorting and degradation; and amino acid metabolism.Table 3Unigenes annotated to the seven databasesComponentNumber of GenesPercentage (%)Annotated in NR52,22962.03Annotated in NT41,03448.73Annotated in KO21,23425.21Annotated in SwissProt41,87749.73Annotated in PFAM40,64048.26Annotated in GO41,64449.45Annotated in KOG15,64518.58Annotated in all Databases909310.79Annotated in at least one Database56,96267.65Total Unigenes84,198100Fig. 2Unigenes enriched in the hierarchical levels of GO, KOG and KEGG terms. The top 20 terms of each classification are included with their corresponding numbers of unigenes

Identification and functional analysis of DEGs {#Sec6}
----------------------------------------------

The read counts of each unigene were sorted by RSEM and transformed to the expected number of fragments per kilobase of transcript sequence per millions base pairs sequenced (FPKM) to characterize the expression levels of the genes. On average, the FPKM value of 73.64% of the unigenes in four libraries was greater than 0.3, which was the minimum value of the subsequent gene expression level analysis. Moreover, the FPKM intervals were 15--60 and \> 60, which signified that the medium and high levels of gene expression were 8.11 and 2.71%, respectively. Using the criterion of padj\< 0.05 for filtering, 1314 DEGs were identified, among which 513 and 801 DEGs were up- and downregulated, respectively (D_R vs L_R, respectively; Fig. [3](#Fig3){ref-type="fig"}). All the DEGs were generated by comparing the expression of unigenes in soil sites D to L (D Vs. L). The up/down regulated DEGs referred to the higher/lower expression of unigenes in plant of soil site D.Fig. 3Analysis of significant DEGs between D_R and L_R. The red and green dots indicate the expression of unigenes that are up- and downregulated in the D_R sample. The volcano plot shows the log~2~(fold change) value, while -log~10~(padj) values are shown on the ordinate. Smaller corrected *p*-values (padj) are correlated with greater corresponding -log~10~(padj) values. Therefore, the points in the upper-left and upper-right sections represent the down- and upregulated genes whose expression levels differed significantly, respectively. With the criterion of padj\< 0.05 for filtering, the red and green dots in the figure are greater than 1.3, and 1314 DEGs were identified. Among them, 513 and 801 DEGs were up- and downregulated, respectively, between D_R and L_R

The DEGs clearly enriched in each GO term are shown in Fig. [4](#Fig4){ref-type="fig"}. DEGs were enriched in 19 BP and MF terms, and the top three enriched terms were metabolic process (BP), catalytic activity (MF), and single-organism process (BP); these terms were associated with 746, 670 and 575 DEGs, respectively, which were equal to percentages of 69.46, 62.38, and 53.54%, respectively. Directed acyclic graphs (DAGs) can intuitively reveal the GO terms of enriched DEGs and their hierarchical relationship. We performed enrichment analysis for the up- and down-regulated DEGs based on their BP and MF, respectively.Fig. 4Illustration of the passion fruit root DEG enrichment of 19 terms with BP and MF categories. In the figure, "oxidoreductase activity, acting on" indicates "oxidoreductase activity, acting on paired donors, with oxidation of a pair of donors resulting in the reduction of molecular oxygen to two molecules of water"

The DAGs of the upregulated DEGs (Fig. [5](#Fig5){ref-type="fig"}) revealed 4 and 8 levels for MF and BP, respectively. Cytokine binding (GO: 0019955) in level 4 was the only term classified as MF. For BP, six terms were enriched in four sublevels. Signal transduction by protein phosphorylation (GO: 0023014) and regulation of shoot system development (GO: 0048831) at level 8 and embryonic morphogenesis (GO: 0048598) at level 6 were enriched separately without enriched parent nodes. Seed germination (GO: 0010029) was a sub-node of both seed germination (GO: 0009845) and regulation of seedling development (GO: 1900140) at level 7; these two terms were subnodes of seedling development (GO: 0090351) at level 6.Fig. 5DAG showing the results of DEG enrichment according to GO analysis. The top ten enrichments were selected in each category of two GO classifications. The main node is described in the boxes, and deeper colors reveal a higher degree of enrichment. Terms with low intermediate expression are not shown. The DAGs were divided into four categories: upregulation in BP (**a**) and MF (**b**) and downregulation in BP (**a**) and MF (**b**)

The downregulated DEGs of 10 terms were enriched in each DAG, which were layered with 7 and 8 levels of MF and BP, respectively. For MF, except level 6, the DEGs were enriched in terms at every level. Three branches were significantly enriched, representing oxidoreductase activity, with the parent node of catalytic activity (GO: 0003824); the second was iron ion binding (GO: 0005506), and the third was tetrapyrrole binding (GO: 0046906) to heme binding (GO: 0020037). With respect to the BP ontology, the DEGs were significantly enriched in the oxidation-reduction process (GO: 0055114); their parent node was associated with the single-organism metabolic process (GO: 0044710). There were also some DEGs enriched in terms involved in cellular amino acid metabolic processes.

Following analysis of the clearly enriched DEGs, the KEGG classification results were presented as a scatter plot (Fig. [6](#Fig6){ref-type="fig"}). For significantly enriched DEGs, the top three pathways were plant hormone signal transduction (31 DEGs), phenylpropanoid biosynthesis (30 DEGs), and biosynthesis of unsaturated fatty acids (26 DEGs), with Q-values of 0.025, 2.78e-06 and 2.81e-09, respectively. The pathways of alpha-linolenic acid metabolism (20 DEGs), phenylalanine metabolism (16 DEGs), ubiquinone and other terpenoid-quinone biosynthesis (11 DEGs), fatty acid degradation (16 DEGs), peroxisomes (20 DEGs), and cysteine and methionine metabolism (20 DEGs) were significantly enriched. The specific information of DEGs enriched in KEGG pathways were showed in supplement file (Additional file [5](#MOESM5){ref-type="media"}: Table S1).Fig. 6The degree of enrichment was measured by the rich factor, Q-value and numbers of enriched genes in the pathways. The top 20 pathways with the highest degree of enrichment are shown. The rich factor is the ratio between the number of enriched DEGs and the number of all annotated genes in the pathway. The Q-value indicates the p-value after multiple hypothesis test correction and ranges from 0 to 1. A higher enrichment is achieved when the Q-value approaches 0

The KEGG enrichment analysis revealed that nine pathways were significantly enriched, with the top three representing plant hormone signal transduction, phenylpropanoid biosynthesis and biosynthesis of unsaturated fatty acids. Thirty-one DEGs were enriched in six sub-pathways of plant hormone signal transduction (with 568 background unigenes): 6 DEGs were enriched in auxin signal transduction, 8 in cytokinin signal transduction, 3 in gibberellin signal transduction, 5 in abscisic acid (ABA) signal transduction, 3 in ethylene (ETH) signal transduction, and 6 in salicylic acid (SA) signal transduction (Fig. [7](#Fig7){ref-type="fig"}). Thirty of the DEGs enriched involved in phenylpropanoid biosynthesis pathway (Additional file [1](#MOESM1){ref-type="media"}: Figure S1) and 26 involved in the biosynthesis of unsaturated fatty acids pathway were downregulated (Additional file [2](#MOESM2){ref-type="media"}: Figure S2).Fig. 7KEGG pathway enrichment analysis of plant hormonal signal transduction pathways. DEGs were enriched in 6 of the 8 subpathways. The red highlighting indicates that a DEG was upregulated, the green indicates that a DEG was downregulated, and yellow indicates both

qPCR validation of DEGs {#Sec7}
-----------------------

We randomly selected 13 DEGs for qPCR assessment to validate the accuracy of the Illumina RNA sequencing (RNA-Seq) results, and DEGs mentioned in discussion were all verified. These DEGs presented an abundance of increased or decreased transcripts between root samples from the two different soils. For qPCR, the level of gene expression was determined by the relative quantity (RQ). Compared to RNA-Seq, the RQ results were consistent with the expression patterns of target genes assessed via qPCR. In addition, we used fold changes to quantify the relevance of the RNA-Seq expression levels and qPCR results. Figure [8](#Fig8){ref-type="fig"} shows a scatter plot of the RNA-Seq data and the qPCR results as log~2~(fold changes), indicating significant similarity (R^2^ = 0.9367, *p* \< 0.01) between the results of the two methods and verifying the accuracy of the RNA-Seq. On the other hand, the trend of RQ results of DEGs listed in discussion in two soil sites was mostly consistent with the prediction of RNA-Seq (Additional file [3](#MOESM3){ref-type="media"}: Figure S3 and Additional file [4](#MOESM4){ref-type="media"}: Figure S4).Fig. 8Scatter plot based on the log~2~(fold changes) of the RNA-Seq and qPCR data. The abscissa represents the log~2~(fold change) value of the read counts from RNA-Seq, while the ordinate indicates the log~2~(fold changes) of the negative delta-delta Ct values from qPCR. The value of R^2^ = 0.9367 (*p* \< 0.01) indicates the significant similarity of the results between RNA-Seq and qPCR, which demonstrated the accuracy of RNA-Seq

Discussion {#Sec8}
==========

DEGs involved in nutrient uptake, transport, and transduction {#Sec9}
-------------------------------------------------------------

Roots play a very important role in plant anchorage, water and nutrient uptake, and the synthesis and storage of metabolites \[[@CR13]\]. N, P, and K are essential nutrient elements and play an indispensable role in plant growth and development. In plants, P is an important component of metabolites and macromolecules, such as ATP, phospholipids and nucleic acids, and is involved in many biological pathways, such as those related to gene expression and signal transduction \[[@CR14]\]. Orthophosphoric acid (H~2~PO~4~^−^, HPO~4~^2−^ and PO~4~^3−^) is the main form of absorption in plants, and there is a severe shortage of P in agricultural soils worldwide. The P starvation response aims to achieve the self-balance of P via a series of physiological and biochemical activities, including the absorption, transportation and reuse of P. Under low P conditions, plant roots activate insoluble P in the soil by secreting organic acids (OAs) to release Pi from refractory P compounds in the soil and increase the content of AP \[[@CR15]\]. Previous studies have shown that the activity of phosphoenolpyruvate carboxylase (PEPC) was increased under conditions of P deficiency, providing a carbon source for the synthesis of OAs \[[@CR16], [@CR17]\]. The results of those studies showed that 27 unigenes were homologous with the PEPC family protein, and one DEG of phosphoenolpyruvate carboxylase kinase 2 (PPCK2) was significantly upregulated. The activity of PEPC is regulated mainly by phosphoenolpyruvate carboxylase kinases (PPCKs) \[[@CR18]\]. The low phosphatase conditions in the soil of the L plots probably stimulated the metabolic activities in the passion fruit plants and promoted the synthesis of OAs for the transformation of insoluble P in the soil.

Adaptation to low P environments and P procurement from the environment require a complex signaling network in plants \[[@CR19], [@CR20]\], which is dependent on the P transport system. P transporters (PHT) play an important role in the absorption and transportation of P \[[@CR21]\] and are grouped into five families: PHT1, PHT2, PHT3, PHO1 and PHO2. In this study, 137 unigenes, associated to 28 PHT1, 3 PHT2, and 25 PHO1, were homologous to those of P transporter families. Among those unigenes, the expression of PHO84, a DEG encoding the protein of a high-affinity P transporter \[[@CR22]\], was significantly higher in the roots from the L site soil than in the roots from the D site soil (Table [4](#Tab4){ref-type="table"}). The results from the present study showed that the amount of AP in the L site soil (0.0026 g/kg) was considerably lower than that in the D site soil (0.2115 g/kg). However, the P content in plants grown at the two sites was approximately equal. The expression of P transporters was inconsistent with the trend of element changes in both the plants and soils. The upregulation of P transporter genes improved the absorption of P by roots at the site L soil.Table 4The characteristics and expression of DEGs associated with nutrient elements uptake, utilization, and signal regulationGene_idAnnotationFPKM(D_R)FPKM(L_R)D_R vs L_RCluster-7898.36820*PHO84*21.89557.785DownCluster-7898.4433*APS1*3.8915.175DownCluster-7898.45233*PAP1*4.68513.32DownCluster-7898.18934*PAP1*0.2253.12DownCluster-7898.27691*PPCK2*12.3239.205DownCluster-7898.61262*POT5*6.9833.955DownCluster-7898.5070*AKT1*2.7759.135DownCluster-7898.5365*HAK5*0.7556.74DownCluster-7898.26747*CLC-b*11.792.735UpCluster-7898.53478*CLC-b*3.5350.355UpCluster-7898.46597*NIPA*14.0148.875DownCluster-7898.49803*CML11*78.065414.37DownCluster-7898.61829*CBL7*2.35519.595DownCluster-7898.46341*GLR2.7*0.543.12DownCluster-7898.5855*GRXC6*2.7627.335DownCluster-7898.61518*CIPK14*1.5910.925DownCluster-7898.2597*CIPK8*4.44255.4625DownThe functional information of DEGs listed were showed in supplemental files (Additional file [6](#MOESM6){ref-type="media"}: Table S2)

In addition, Pi starvation induces the synthesis of acid phosphatases (APases) in plants. In roots, APases are important for the hydrolysis and mobilization of Pi from extracellular phosphomonoesters for plant nutrition \[[@CR23]\]. In the present study, 116 unigenes were annotated to APase genes, and one *ASP* and two *PAP1* genes were significantly upregulated in the roots from the site L soil (Table [4](#Tab4){ref-type="table"}). A low P environment reportedly induces the transcription of the *PAPs* \[[@CR24]\]. The results of the present study showed that *P. edulis* could increase the activity of APases to improve the efficiency of P utilization and the ability of P absorption to adapt to a low P stress environment. These responses may play an important role in the stability of P balance in plants.

Plants contain a large amount of K, which exists in an ionic form within cell sap or is adsorbed on the surface of protoplast colloids. K^+^ plays an important role in maintaining the balance of cell potential and osmotic potential, and K is absorbed through K channels and K^+^ transporter in plants \[[@CR25]\]. In this study, 219 unigenes were annotated to K transport-related proteins, including 111 K channel genes and 79 K^+^ transporter genes; among the K^+^ transporter unigenes, 15 were annotated to high-affinity K^+^ transporters (HAKs). The amount of K in the root samples at site D was clearly greater than that at site L. However, DEG analysis revealed that the expression of the *HAK5*, *POT5*, and *AKT1* genes was downregulated. The expression model was not consistent with the changes in K in the plants and soils at sites D and L (Table [4](#Tab4){ref-type="table"}). Previous studies showed that *HAK5* \[[@CR26]\] and *AKT1* \[[@CR27]\] were inducible in a high-affinity K^+^ uptake system under low K conditions, which may suggest that the low concentration of K in the site L soil mediates the expression of these genes for K uptake.

Alkali hydrolyzable N consisted of inorganic N (ammonium \[NH~4~^+^\]-N, nitrate \[NO~3~^−^\]-N) and hydrolytic organic N (amino acids, amides and hydrolyzable proteins). NH~4~^+^-N and NO~3~^−^-N are the two major N sources for higher plants, and NO~3~^−^ is the most abundant source of N \[[@CR28]\]. NH~4~^+^ AMT transporters are the major pathways for the uptake of NH~4~^+^ in roots \[[@CR29]\]. In the present study, 29 unigenes were homologous to NH~4~^+^ transporters, including 16 AMT type genes. However, the expression of these unigenes was not markedly different between the roots at the two soil sites. On the other hand, four gene family members were shown to be involved in NO~3~^−^ uptake, allocation, and storage in higher plants: NO~3~^−^ transporter 1/peptide transporter (NRT1/PTR), NRT2, chloride channel (CLC), and slow anion channel-associated 1 homolog 3 (SLAC1/SLAH) \[[@CR30]\]. With respect to the unigene annotations, 75 unigenes were homologous to NO~3~^−^ transport-related genes, including 61 *NRT1*/*PTR*s, 4 *NRT2*s, and *5 CLC*s. Differential expression analysis revealed that 2 *CLC-b* genes were expressed at higher levels in roots from site D than in roots from site L (Table [4](#Tab4){ref-type="table"}). In addition, the results showed a tendency toward similarity in the amount of AN in the soil; however, the concentration differed between plants, with a clear increase in plant growth at site D compared with that at site L. The uptake of NO~3~^−^ in plants is influenced by many factors; NO~3~^−^-N is absorbed faster in lower-pH soil solution, and nutrient elements, such as P \[[@CR31]\] and K \[[@CR32]\], promote the absorption of N. The results possibly indicated that the relatively low pH and relatively high concentrations of P, K, and Ca may have increased the expression of NO~3~^−^ transport-related genes in the roots at site D, influencing the intracellular exchange of NO~3~^−^.

Mg^2+^ plays an important role in the physiological processes of plant growth and development. For plant adaptation to different soil conditions, members of two Mg transporter families have been reported to regulate intra- and extracellular Mg^2+^ transport in higher plants: CorA-like Mg^2+^ transport systems \[[@CR33]\] and Mg/proton exchangers \[[@CR34]\]. In the present study, 93 unigenes were annotated to Mg^2+^ transport, among which 26 CorA-like Mg^2+^ transporters and 5 Mg/proton exchangers were identified. One DEG encoding the Mg^2+^ transporter *NIPA* was downregulated in the roots of the soil at site D compared with that at site L (Table [4](#Tab4){ref-type="table"}). Moreover, the concentration of Mg^2+^ in the plants and soils was higher at site L than at site D, and the expression pattern was in accordance with the trend of the Mg^2+^ content in the roots and soils, which may indicate that the gene expression level tended to be positively regulated by the changes in elements.

As a universal signaling messenger, Ca^2+^ is an essential nutrient element of plants and regulates many physiological processes in plant cells \[[@CR35]\]. Hyperpolarization-activated Ca^2+^ channels (HACCs), depolarization-activated Ca^2+^ channels (DACCs,), Ca^2+^-permeable outward-rectifying K^+^ channels (KORCs), and voltage-insensitive cation channels (VICCs) are four important Ca^2+^ channels in the plasma membrane \[[@CR36]\]. In the present study, 134 unigenes were annotated to Ca^2+^ channel genes, including 30 annexin, 2 *TPC1*, 6 *SKOR*, 1 *GORK*, 41 *CNGC*, and 52 *GLR* genes. Of these genes, those coding for Ca^2+^ channels were highly expressed under high concentrations of Ca ions in the soils at both sites. Specifically, 59 unigenes were filtered after they were standardized to FPKM\> 3; among them, 6 unigenes with an FPKM value \> 60 were significantly expressed, and 27 unigenes were filtered by 15 \< FPKM\< 60. The highly expressed unigenes related to Ca^2+^ channels suggested that the concentration of Ca^2+^ in the cytoplasm was increased under high Ca concentrations.

The soils of Karst areas are based on dissolved Ca-rich parent rock, and the native plants in those regions have evolved a special Ca ion adaptation system \[[@CR4]\]. Therefore, this system could be the source of the adaptation mode for 'Pingtang 1' plants growing on Karst landforms. Roots take up Ca via Ca^2+^ channels and transport it through Ca^2+^-ATPase and Ca^2+^/H^+^ antiporters to maintain low concentrations of Ca in the cytoplasm \[[@CR30]\]. Plants that have adapted to Ca-rich environments by enriching the amount of Ca ions can transport those ions and store them via Ca^2+^/H^+^ reverse transporters and Ca^2+^-ATPases. Our results showed that 113 unigenes were Ca^2+^-ATPase genes; among these genes, 33 were filtered according to the FPKM\> 3 standard, while 8 were highly expressed, with FPKM values\> 15. Moreover, 26 unigenes were homologous to Ca^2+^/H^+^ antiporters, and 9 were highly expressed, with FPKM values \> 3. These findings suggest that the high expression levels of Ca^2+^-ATPase and Ca^2+^/H^+^ antiporters may be responsible for maintaining low Ca concentrations in the cytoplasm.

In response to varied extracellular Ca^2+^ levels, relatively constant cellular concentrations of Ca^2+^ suggested a regulatory mechanism for Ca^2+^ homeostasis \[[@CR37]\]. Ca^2+^-sensing receptors (CASs), which are localized at the plasma membrane, may be primary transducers of extracellular Ca^2+^ in plants \[[@CR38]\]. In this study, 7 unigenes were homologous to CASs, and 3 of these genes, whose FPKM values were \> 3, were expressed. Ca-dependent protein kinase (CDPK), calmodulin (CaM) and CaM B-like protein (CBL) are three major Ca^2+^ receptors involved in signal transduction. CDPKs comprise four types of subfamilies: CDPKs, CDPK-related proteins (CRKs), CaM-dependent protein kinases (CaMKs), and Ca^2+^ and CaM-dependent protein kinases (CCaMKs). CaM transduces Ca^2+^ signals through CaM-binding protein (CaMBP) \[[@CR39]\], and CBL transduces Ca^2+^ signals only by interacting with CBL-interacting proteins (CIPKs) \[[@CR40]\]. Our results revealed 106 predicted CaMs, 32 CaMKs, 67 CDPKs, 29 CDMKs, 16 CRKs, 11 CCaMKs, 24 CBLs and 51 CaMBPs in the libraries. Of these genes, 67 were filtered according to the FPKM\> 3 standard. The expression patterns of the Ca^2+^ receptor genes revealed that high Ca^2+^ concentrations in the cytoplasm mediated downstream reactions in cells. On the other hand, the expression of DEGs annotated to Ca^2+^ channels (*GLR2.7*, *GRXC6*), Ca^2+^ signal transport proteins and downstream receptors (*CIPK14*, *CIPK8*, *CML11*, *CBL7*) was downregulated in the roots of the soil at site D compared with that at site L (Table [4](#Tab4){ref-type="table"}). Considering the consistency of Ca^2+^ changes in plants and soils, Ca^2+^ absorption and utilization in cells may tend to be positively regulated by elemental changes in the soil environment.

In this study, at the transcriptional level, differential analysis of element transporter genes revealed that the uptake, transport and utilization of nutrient elements were affected by the elemental conditions at the soil sites. The changes in the amounts of elements in the plants were not consistent with the changes in the soils at the two sites, which may indicate that interactions occurred between elements when the plants responded to the different soil environments. Furthermore, the expression patterns of DEGs involved in the plant hormone signal transduction pathway (Fid. 7) suggested that the responses of plants were different in the soils of different nutrient element contents. The expression pattern of genes related to Ca^2+^ transport and transduction may suggest that Ca-rich conditions induce a series of responses concerning Ca uptake, transport, and downstream reactions of Ca^2+^ receptors, which occur to maintain a dynamic Ca balance within 'Pingtang 1' plant cells to allow adaption to Karst areas. These results will contribute to revealing the adaptive mechanism.

DEGs enriched in top-three significant KEGG pathways {#Sec10}
----------------------------------------------------

In the study, 'Pingtang 1' was a specifically cold-resistant variety of *P. edulis* in karst areas. It was necessary to research the mechanism of adaptability to low temperature in two soil sites of karst areas. From the results of DEGs in KEGG enrichment, we focused the sight on the top-three significant pathways, which were plant hormone signal transduction, phenylpropanoid biosynthesis, and biosynthesis of unsaturated fatty acids. On this basis, combining with the mineral elements condition of soil and plant in two soil sites, there were some analysis and suggestions.

In previous researches, cold-resistances of plants were associated with the components of the fatty acid of the plasma membrane, the stability of cell structure would be improved by increasing membrane lipid unsaturation \[[@CR41]\]. DEGs in the pathway of biosynthesis of unsaturated fatty acids were all down-regulated. Several DEGs were homologous with FAD2, which encodes the enzyme that is essential for polyunsaturated lipid synthesis in the endoplasmic reticulum \[[@CR42]\]. Hence, the differential expression of these DEGs related to lipid unsaturation may influence the responses to low temperature of plants.

As for phenylpropanoid biosynthesis pathway, DEGs involved in the biosynthesis of cinnamaldehyde, p-coumaryl alcohol, caffeylalcohol, 5-Hydroxyconiferyl alcohol, and sinapyl alcohol were obviously down-regulated, however, DEGs encoded peroxidase \[EC:1.11.1.7\] were both up-regulated and down-regulated. Those alcohols were involved in the biosynthesis of secondary metabolism, the secondary metabolic process was believed to be the result of adaptation in dealing with the relationship between plants and the ecological environment. And therefore, the down-regulated DEGs may decrease the synthesis of related secondary products and further affect potential adaptability of *P. edulis* in karst areas. Research has shown exogenous K elements would diminish the expression of gene 4CL and F5H \[[@CR43]\]. And the amount of AK of soil and tissues in site D were lower than D, which may indicate the association with significantly down-regulated expression of 4CL and F5H.

In the pathway of plant hormone signal transduction, auxin regulates plant growth and development by altering the expression of diverse genes \[[@CR44]\], while DEGs related to three auxin-responsive AUX/IAA and SAUR families were both up-regulated, and SAUR and GH3 were down-regulated in the results. In cytokinin signaling, the type-B response regulators mediate cytokinin signal transduction (B-ARR) \[[@CR45]\] were down-regulated. And the other type response regulator (A-ARR) was up-regulated, which was reported that could negatively regulate cold resistance of plants \[[@CR46]\]. In the research, DEGs of ARR-A would possibly induce the further responses of chilling-resistance of 'Pingtang 1' in soil site L. In the downstream pathway of Gibberellin (GA), DELLA, as a receptor of GA \[[@CR47]\], the higher expression of which in soil site D may suggest the positive impact to regulatory effects of GA in site D's plant. Abscisic acid inhibits type 2C protein phosphatases (PP2C) via the PYR/PYL family, and high PP2C activity may prevent phosphorylation and activation of SnRK2s and downstream factors \[[@CR48]\]. And moreover, SnRK2s could phosphorylate the ABA-responsive element binding factor family to activate ABA-responsive genes \[[@CR49]\], which would further encode protective proteins and enzymes for osmolyte synthesis for stress signal transduction in plants \[[@CR50]\]. DEGs related to PYL and SnRK2 were down-regulated, while PP2C were up-regulated in abscisic acid sub-pathway. The biosynthesis of ethylene was induced by abiotic stress, such as chemicals, metals, and extreme temperatures \[[@CR51], [@CR52]\]. In signal pathway leading from salicylic acid (SA), NPR1 was a regulatory linkage between TGA factors and PR genes \[[@CR53]\]. Up-regulated DEGs annotated to NPR1 and RP-1 proteins were and down-regulated TGA may suggest us the differential responses to stress for plant defense in two kinds of soil sites. Comprehensively, pathways of plant hormone signal transduction, phenylpropanoid biosynthesis, and biosynthesis of unsaturated fatty acids participate in complicated physiological processes of plants. The expression pattern of these DEGs may influence the downstream responses about plant growth, development, adaptivity to the environment, and the resistance to stress.

'Pingtang 1' were reported both adaptive to the cold condition and different soil types in karst areas in southwestern of China \[[@CR10]\], in our study, the physiological characteristics of two cutting plants grew in two soil sites were similar after carefully observation. Many factors existed in microenvironments would influence the plant growth and development, and the effect from nutrient elements on the adaptation of *P. edulis* to environmental stress was partly limited. In view of the impact of the environmental condition such as terrain, light, annual precipitation, and human activity, the two adjacent soil sites with different soil types were selected for minimizing the interference of other factors while designing the experiment. On this basis, combined with the connection of different elements concentration of soil and plant in two site plots, the key point was concentrated on the analysis of molecular mechanism about the adaptability of passion fruits in karst areas. However, the influence of other environmental factors was inevitable, this research was our first step to study the adaptive mechanism of passion fruits in karst areas, and we would like to consider more factors in the next researches for further exploration.

Conclusion {#Sec11}
==========

'Pingtang 1', a new variety of *P. edulis*, which is highly adaptable in Karst areas. After integrated analysis upon nutrient elements differences in both the soils and plants and transcriptomic profile in plant roots, the results of DEGs expression patterns related to utilization, transport and transduction of mineral elements, and further analysis about DEGs in top-three enriched pathways will support for understanding to exploring the adaptive mechanism of plants.

Methods {#Sec12}
=======

Plant materials and growth conditions {#Sec13}
-------------------------------------

'Pingtang 1', the cold-tolerant variety of *P. edulis*, was massively propagated by cuttings and served as the research material. Cuttings of 'Pingtang 1' were cultivated in L rocky desertification areas and sandy D rocky desertification areas. The two plots are located in Kedu town, Pingtang Country, Guizhou Province (25°43′0.53″N, 106°48′19.48″E, altitude 852 m above sea level, mid-subtropical humid monsoon climate). Soil samples were collected from three randomly selected locations in each test plot. The roots, stems and leaves of 2-year-old plants were sampled as three biological replicates per test plot. Two root samples selected from those three biological replicates for RNA extraction were immediately frozen in liquid nitrogen (N) and stored at − 80 °C in an ultralow temperature freezer. We declare plant materials (*P. edulis*) and the field culture used in the research complied with government regulations.

Determination of nutrient elements {#Sec14}
----------------------------------

Soil nutrient elements were detected as described in the most recent Agriculture Industry Standard of China (2015). Briefly, the soil pH was measured in a water:soil solution (1:2.5) via a pH meter. AP was extracted in accordance with the molybdenum antimony colorimetry method at a soil solution ratio of 1:10, and AN was determined using the alkaline hydrolysis diffusion method after reduction with FeSO~4~ and Zn. With respect to the medicinal preparation of ammonium acetate, AK was extracted by the flame photometry method, while ECa and EMg were extracted by the atomic absorption spectrophotometry method.

The proportions of the N, K, P, Ca, and Mg elements in the roots, stems, and leaves were determined in 2-year-old plants growing in the two plots. After H~2~SO~4~-H~2~O~2~ digestion, the N content was determined by the semimicro alkali distillation method, while the molybdate yellow colorimetric method and the flame spectrophotometric method were used to determine the contents of P and K, respectively. The atomic absorption spectrophotometry method was used to determine the total Ca, and Mg was determined after digestion by dry ashing. The whole-plant proportion of each element was calculated by the weighted average of three sets of samples.

RNA isolation, cDNA library construction and sequencing {#Sec15}
-------------------------------------------------------

The total RNA of the roots was isolated using a RNeasy Plant Mini Kit (Qiagen, Hilden, Germany). The RNA quality was monitored via 1% agarose gel electrophoresis and via an Agilent Bioanalyzer 2100 (Agilent Technologies, CA, USA). The cDNA libraries of the L plot were named L_R, and those of the D plot were named D_R.

The four cDNA libraries were generated using an NEBNext Ultra RNA Library Prep Kit (NEB, USA) according to the manufacturer's instructions (Illumina). According to the manufacturer's recommendations, clustering of the index-coded samples was performed on a cBot Cluster Generation System using TruSeq PE Cluster Kit v3-cBot-HS (Illumina, USA). The prepared libraries were sequenced, and 150-bp paired-end reads were generated (Illumina HiSeq 2500, USA). The raw images were transformed using CASAVA base-calling in the FASTQ format (raw data), and all raw reads were deposited into the National Center for Biotechnology Information (NCBI) Sequence Read Archive (SRA accession no. SRP150688).

Transcriptome assembly and annotation {#Sec16}
-------------------------------------

Clean data (clean reads) were generated by filtering adapter sequences and low-quality reads that consisted of more than 10% poly-N or contained more than 50% of bases that had a Q-value ≤20. All downstream analyses were based on high-quality clean data. Transcriptome assembly was achieved using Trinity version r20140413p1 \[[@CR54]\] based on left.fq and right.fq, with the min_kmer_cov value set to 2 by default and all other parameters set to their default. The clean reads of various isoforms of one gene were assembled into distinct transcripts, but the same subcomponent, which could be considered a gene, and the longest transcript of each subcomponent were defined as "unigenes" for annotation. All assembled unigenes were queried via BLAST against the Nr, Nt, Pfam, Clusters of Orthologous Groups of proteins (KOG/COG), SwissProt, Kyoto Encyclopedia of Genes and Genomes (KEGG) Orthology (KO) and Gene Ontology (GO) databases using the Blast2GO program with a cutoff E-value of 10^− 6^.

Differential expression analysis {#Sec17}
--------------------------------

The relative expression levels of the unigenes were estimated by mapping clean reads to the Trinity transcriptome assembly using RNASeq by Expectation Maximization (RSEM) \[[@CR55]\] with the trimmed mean of M-values method. Differential expression analysis of the unigenes was performed using read counts with the DESeq R package. The resulting *p*-values were adjusted to Q-values to compensate for multiple hypothesis testing. Unigenes with an adjusted Q-value \< 0.05 (reported by DESeq) were assigned as DEGs. A volcano plot of the DEGs was constructed using the log~10~(padj) and log~2~(fold change) values. GO and KEGG enrichment analyses of the DEGs were implemented using the GOseq R package and KOBAS software \[[@CR39]\], respectively.

Quantitative real-time polymerase chain reaction (qPCR) {#Sec18}
-------------------------------------------------------

Total RNA was isolated from various samples and treated with RNAse-free DNase I. The DNase-treated RNA was reverse transcribed using random hexamer primers; the specific primers were designed to generate amplification fragments of 70--150 bp. Quantitative real-time polymerase chain reaction (qPCR) was performed on an ABI ViiA 7 Real-Time PCR platform using FastStart Universal SYBR Green Master with ROX in accordance with the manufacturer's protocol; all reactions were performed in triplicate. The relative expression levels were calculated using the delta-delta Ct method, and *HIS* was used as a housekeeping gene \[[@CR10]\].

Additional files
================

 {#Sec19}

Additional file 1:**Figure S1.** Down-regulated DEGs enriched in phenylpropanoid biosynthesis pathway. (JPG 343 kb) Additional file 2:**Figure S2.** Down-regulated DEGs involved in the biosynthesis of unsaturated fatty acids pathway. (JPG 220 kb) Additional file 3:**Figure S3.** The validation of DEGs listed in Table [4](#Tab4){ref-type="table"}. The figure combined bar graph of gene expression and the scatter diagram of correlation analysis between RNA-seq and qPCR. DEGs involved were mentioned in Table [4](#Tab4){ref-type="table"}. (JPG 469 kb) Additional file 4:**Figure S4.** The validation of DEGs in KEGG pathways. The figure combined bar graph of gene expression and the scatter diagram of correlation analysis between RNA-seq and qPCR. DEGs involved were enriched in the pathways of biosynthesis of unsaturated fatty acids, phenylpropanoid biosynthesis and plant hormone signal transduction, respectively. (JPG 392 kb) Additional file 5:**Table S1.** The annotation of DEGs listed in Table [4](#Tab4){ref-type="table"}. (XLS 55 kb) Additional file 6:**Table S2.** The details of DEGs enriched in KEGG pathways. (XLS 28 kb)

4CL

:   4-coumarate-CoA ligase

AK

:   Available potassium

AN

:   Available nitrogen

AP

:   Available phosphorus

CAD

:   Cinnamyl-alcohol dehydrogenase

DAGs

:   Directed acyclic graphs

DEGs

:   Differentially expressed genes

ECa

:   Exchangeable calcium

EMg

:   Exchangeable magnesium

F5H

:   Ferulate-5-hydroxylase

FPKM

:   Fragments per kilobase of transcript sequence per millions base pairs sequenced

NGS

:   Next-generation sequencing

OAs

:   Organic acids

I would like to thank Dr. Ying Yang for the suggestions when analyzing the data, and appreciate for the help of professor Yousry from UBC during the writing.

Funding {#FPar1}
=======

This project was supported by the following programs: the design of the study was supported by all the programs; the collection was supported by the Key Science and Technology Program of Guizhou Province ((2016)2525), the National Key Research and Development Program of China (2016YFC050260404), the Key Science and Technology Program of Guangzhou Province (201604020006), the Key Science and Technology Program of Guizhou Province ((2017)5720--003), the Key Science and Technology Program of Guizhou Province ((2017)2570--2) and the Natural Science Foundation of Guizhou Province (\[2018\]1148); the analysis, interpretation of data, and the writing was supported by the Priority Academic Program Development of Jiangsu Higher Education Institutions (PAPD) and the Collaborative Innovation Plan of Jiangsu Higher Education (CIP).

Availability of data and materials {#FPar2}
==================================

The data of RNA-Seq is available in NCBI database (SRA accession no. SRP150688). The datasets generated and analyzed of the current research are available from the corresponding author on reasonable request.

MX conceived and designed the experiments. ADL involved in the research design. YT collected the samples. MXX executed the experiments. ZMS helped the experiments. MXX analyzed the results and wrote the paper. All the authors have read and approved the final manuscript.

Ethics approval and consent to participate {#FPar3}
==========================================

Not applicable.

Consent for publication {#FPar4}
=======================

Not applicable.

Competing interests {#FPar5}
===================

The authors declare that they have no competing interests.

Publisher's Note {#FPar6}
================

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
