Reference Gene Selection for RT-qPCR Assays in Different Tissues of Huperzia Serrata Based on Full-Length Transcriptome Sequencing

Huperzia serrata (H. serrata) produces various types of effectively lycopodium alkaloids, especially Huperzine A (HupA), which is a promising drug for the treatment of Alzheimer’s disease. Numerous studies focused on the chemistry, bioactivities, toxicology, and clinical trials of HupA, however, the public genomic and transcriptomic resources are very limited for H. serrata research, especially for the selection of optimum reference genes. Based on the full-length transcriptome datasets and previous studies, thirteen candidate reference genes were selected in different tissue of H. serrata. Then, two optimal reference genes GAPDHB and HisH2A were conrmed by various analysis softwares. In order to further verify the accuracy of the two reference genes, they were used to analyze the expression patterns of four HupA-biosynthetic related genes (lysine decarboxylas, RS-norcoclaurine 6-O-methyltransferase, cytochrome P45072A1 and copper amine oxidase). The data suggested that the expression trends of HupA-biosynthetic related genes were consistent with them in transcriptome sequencing in different tissue of H. serrata. obtained between transcriptome sequencing and RT-qPCR when researched the expression prole of four HupA-biosynthetic genes, lysine decarboxylase (LDC), RS-norcoclaurine 6-O-methyltransferase (MET), cytochrome P45072A1 (CYP) and copper amine oxidase (CAO). This further veried that GAPDHB and HisH2A were suitable for HupA-biosynthetic gene expression normalization. This work provides suitable RGs for the subsequent research of HupA-biosynthetic in H. was the least expressed gene on account of its highest mean Ct value (29.43). All candidate genes showed expression variability in different samples as evident from a wide range of Ct values. Genes such as GAPDHB and EFTS showed relatively smaller variation (< 2 cycles), while others like UBQ11 had the highest expression variation (3.07 cycles). The results indicated that there was still variable expression even for relative stable housekeeping genes.


Background
Huperzia serrata (H. serrata) belongs to the Huperzia genus, Lycopodiaceae order. The whole plant of H. serrata has been used as a medicine in China to treat different kinds of ailments, including bruises, strains, swelling, rheumatism, schizophrenia, myasthenia gravis, and fever since 739 (during the Tang Dynasty) [10]. H. serrata has been widely known as a medicinal plant since Chinese scientists discovered Huperzine A (HupA) from it during the 1980s [15]. HupA is a promising candidate drug for treating Alzheimer's disease (AD), it could improve cognitive function, daily living activity, and global clinical assessment in patients with AD disease, with relatively few and mild adverse effects [34,25]. However, H. serrata is scarce in nature and grows very slowly in specialized habitats. Furthermore, the HupA content is very low in H. serrata [20]. At present, the rapidly growing demand has put H. serrata resources on the brink of extinction. Although, a lot of efforts have been focused on arti cial culture and tissues culture for H. serrata production, the results were unsatisfactory. Now, researchers try to improve HupA content by studying the gene information of HupA biosynthesis. However, the public genomic and transcriptomic resources are very limited. Only two papers focused on transcriptomic resources [2,36]. Real-time quantitative PCR (RT-qPCR) has been widely used in gene expression measurement in transcriptional level. Identi cation of suitable reference genes (RGs) is pre-requisite for RT-qPCR assays [29,3]. Many housekeeping genes have been used as RGs under different experimental conditions, such as actin, tubulin, elongation factor (EF), 18S ribosomal RNA (18S rRNA), glyceraldehyde-3-phosphate dehydrogenase (GAPDH), histone, ubiquitin and so on [21,9]. However, there is no RG suitable to all biological systems because all previous studies have revealed alterations for housekeeping genes expression. So, researchers reach a consensus that speci c RG for a given species and treatment needs to be identi ed rstly [3]. Unfortunately, the reportorial RGs are not suitable for the research of HupA-biosynthetic genes [35].
In present study, in order to obtain the optimal RGs for researching HupA-biosynthetic, we detected the concentration of HupA and carried out full-length transcriptome sequencing for the different tissues of H. serrata. Based on full-length transcriptome sequencing data and previous studies, thirteen candidate RGs were selected. Finally, GAPDHB and HisH2A stood out among the thirteen candidate RGs, and became the best combination for normalization in different tissues of H. serrata by the comprehensive analysis of four softwares. Using GAPDHB and HisH2A as RGs, the consistent expression trend was obtained between transcriptome sequencing and RT-qPCR when researched the expression pro le of four HupA-biosynthetic genes, lysine decarboxylase (LDC), RS-norcoclaurine 6-O-methyltransferase (MET), cytochrome P45072A1 (CYP) and copper amine oxidase (CAO). This further veri ed that GAPDHB and HisH2A were suitable for HupA-biosynthetic gene expression normalization. This work provides suitable RGs for the subsequent research of HupA-biosynthetic in H. serrata.

Plant Materials
H. serrata plants were collected from Hanzhong, Shaanxi, China (107°09′/32°30′) in March 2018. All materials used in this study were identi ed by phytotaxonomist. The plants were rinsed carefully by running water. Root, stem, and leaves were collected in liquid nitrogen and were immediately frozen at -80°C for RNA extraction. Other plant materials were dried at 60°C and powdered for determining HupA content.

HPLC parameters and conditions
HupA were extracted from the different plant tissues as previously described [19,12]. After the plant material was dried and milled, 100 mg each of powdered plant tissues was extracted by adding 2% (2:100, w/v) aqueous tartaric acid (5 ml) for overnight and then sonicating for 2 h at 25°C. Centrifuging for 30 min at RT and the upper extraction solution was ltered into a 1.5 mL measuring ask through a 0.45 μm lter. Finally, the ltered solutions (10 µL) were injected into the HPLC system (LC-20AT, Shimadzu, Japan) for detection of HupA content.

Candidate RGs selection and primers design
RGs for this study were selected from full-length transcriptome sequencing dates and previous studies where these were found to have stable expression in other plants by using local NCBI-blast (version 2.4.0+). Thirteen candidate RGs (MH560040 -MH560049, MZ042627-MZ042629) were selected (Table 1). Gene speci c primers for each RG were designed using the Primer 5.0. Conserved domains of RGs were evaluated and the primer binding positions were presented. Initially, primer speci city was veri ed by RT-qPCR and con rmed with 2% (w/v) agarose gel electrophoresis and melting curve.

RNA extraction and cDNA synthesis
The total RNA was extracted according to the modi ed CTAB method [11]. RNA samples were treated with DNase I (Ambion, Waltham, MA, USA) to remove any DNA contamination. Using cDNA synthesis kit (Roche, Basel, Switzerland), rst strand cDNA was prepared with 3 μg RNA as manufacturer's instructions.

RT-qPCR analysis
The RT-qPCR reactions were performed with FastStart Universial SYBR GreenMaster (Roche, Basel, Switzerland) on a CFX-96 thermocycling system (Bio-Rad, Hercules, CA, USA). Each RT-qPCR reaction was performed as described previously [24]. PCR ampli cations were carried out with the following cycling conditions: one cycle at 95°C (180 s), followed by 40 cycles of denaturation at 95°C (30 s), annealing at 58°C for 10 s and extension at 72°C for 20 s [16]. Finally, melt curve analyses were done by slowly heating the PCR mixtures from 58 to 95°C. Ampli cation e ciencies (E) and correlation coe cients (R 2 ) for each primer pair were calculated by LinRegPCR program [26]. In the negative control group, RT-qPCR was performed using water instead of cDNA as the template. Three technical replicates were analyzed for each biological sample, and each experiment comprised three independent biological replicates.
For geNorm and NormFinder, the raw Ct values were converted into the relative quantities using the formula 2 -ΔCt (ΔCt = each corresponding Ct value -lowest Ct value). M value, which respents average expression stability, was calculated in geNorm. The candidate RGs showing a higher M value (M>1.5) are not considered for normalization studies [29]. geNorm software was also used to con rm the numbers of RGs with pair wise variation (Vn/Vn+1, n refers to the RGs number) [29]. NormFinder provides the stability value for each gene, which is a direct measure of the estimated expression variation enabling evaluation of the systematic error introduced when using the gene for normalization [1]. For BestKeeper, the raw Ct values and ampli cation e ciencies were used to calculate the coe cient of variation (CV) and standard deviation (SD). The most stable genes are determined to be those which exhibit the lowest CV and SD (CV ± SD). The comprehensive ranking order was recommended on the basis of geometric mean (GM) by RefFinder [37].

Validation of RGs
The primer of four Hup A-biosynthetic genes LDC (GO914645) MET (GO914756), CYP (GO914428) and CAO (JN247732) [17,33,28] were designed using the Primer 5.0. The combination of the top two best ranked RGs and worst ranked RGs were used to standardize the expression of two target genes. The target gene expression data was normalized using the geometric mean values calculated for the RG pairs [29]. Relative expression level and fold change were determined using the comparative 2 -ΔΔCt method [22]. One-way analysis of variance was performed using SPSS software (Version 16.0, SPSS Inc., Chicago, IL, USA).

Results And Discussion
The Hup A content analyses HPLC-UV was performed to detect the HupA content in H. serrata. Typical chromatograms from HupA standard and three tested samples are shown in Figure S1, indicating that HupA has good peak shape and is well separated from different tissues. A linear relationship exists between the peak area (measured at 308 nm) and the concentration of HupA in the sample injected into the HPLC. The results showed that there was obvious difference for the HupA content in different tissues. The highest HupA content (72 μg/g) was found in the leaves of H. serrata. The lowest content (19 μg/g) of HupA was found in root tissues of H. serrata ( Figure 1).

The screening of candidate RGs
In consideration of the varied difference of Hup A concentrations in different tissues, the root, stem and leaf samples were collected and proceeded the full-length transcriptome sequencing by Nanopore. After assembly, 43,443 unigenes were retrieved. CPM counts per million is the index for measuring the expression of unigenes. Based on the CPM value and reported literatures, ten traditional RGs and three new RGs candidate were choosed. The three new RGs candidate had stable expression in full-length transcriptome sequencing. They were annotated hypothetical or uncharacterized proteins by NCBI Nr database, furthermore, they were not used as RGs before. The three new RGs candidate as following: ONT.10684 represented the high expression level (CPM over 100), EVM0022608 was the middle level (CPM 29-34), and EVM0017784 was the low level (CPM less than 5). The detail information of total thirteen candidate RGs was showed in Table   1.
Veri cation of the primer speci city and RT-qPCR ampli cation e ciency The primer information of thirteen candidate RGs was given in Table 2. Each primer pair was designed except the conserved domains to ensure the speci city ( Figure S2). Initially, the agarose gel electrophoresis yielded a speci c fragment of expected size ( Figure S3A, S3B and S3C). Further, the melting curve analysis in the RT-qPCR reaction showed the single peak for each primer pair indicating an absence of non-speci c product ampli cation ( Figure S4). For all primer pairs, the ampli cation e ciencies were spanning from 90.4% to 103.6%, and the correlation coe cient (R 2 ) were greater than 0.990 (Table 2). Taken together, these results indicated each primer pair was speci city and the RT-qPCR assays were highly e cient.

Expression pro les of candidate RGs
The expression pro les of RT-qPCR products for all experimental samples are shown in Figure2. The results illustrated that the mean Ct values of all RGs ranged from 24.04 to 29.43. Lower Ct value indicates the higher expression abundance, conversely means the lower expression pro les. EF1dt and UBQ1 were highly expressed with mean Ct values between 24.04 and 24.08 while EFTS was the least expressed gene on account of its highest mean Ct value (29.43). All candidate genes showed expression variability in different samples as evident from a wide range of Ct values. Genes such as GAPDHB and EFTS showed relatively smaller variation (< 2 cycles), while others like UBQ11 had the highest expression variation (3.07 cycles). The results indicated that there was still variable expression even for relative stable housekeeping genes.

geNorm analysis
To identify the most stable RG, geNorm algorithm calculated the average expression stability values (M values) of each RG. As Figure3 shown, each M value was less than 1.5, which suggested the appropriateness of all RGs for normalization consideration in different tissues of H. serrata. Concretely analyzing, EF1dt, HisH2A and GAPDHB were the most stable genes in each H. serrata samples, while HisH3.3 and EFTS were the least stable genes in each H. serrata samples.

NormFinder analysis
NormFinder evaluates each RG according to the stability value. Lower stability value indicates more stable gene expression, and vice versa. As shown in Table 3, GAPDHB and HisH2A were obviously stable in all samples, and EFTS (highest stability value = 0.210) was the least stable gene. For the root samples, HisH2A and a-tub3 were most stable, and the EFTS (stability value = 0.506) still was the least stable gene. Whereas GAPDHB and HisH2A were the most stable gene and HisH3.3 (stability value = 0.358) was the least one in stem. In leaf tissues, the most stable RGs were EF1dt, HisH2A and, GAPDHB, meanwhile, the least stable RG was Actin7 (stability value = 0.364). Overall, with NormFinder analysis, GAPDHB, HisH2A and EF1dt were the most stable genes, while EFTS and HisH3.3 were the least stable genes in different tissues of H. serrata.

BestKeeper analysis
The stability standard deviation (SD) and its relationship to the BestKeeper index were considered as two important evaluation criteria in BestKeeper analysis [23]. The results showed that each RG had a SD value < 1.0, which indicated that the candidate RGs were relatively stable for RT-qPCR normalization. In present, GAPDHB, EF1dt and HisH2A were the top three ranked genes with lowest CV ± SD values in all samples, stem and leaf tissues (Table 4). In root samples, the top three ranked genes were HisH2A, UBQ11 and a-tublin. Rather, HisH3.3 was deemed to the least RG with the highest CV ± SD value (27.35 ±0.21 and 28.66 ±0.20) in all samples and stem tissues, while EFTS (in root tissues) and Actin7 (in leaf samples) showed the least stable expression. Taken together, with BestKeeper analysis, GAPDHB, EF1dt and HisH2A were the most stable genes, while HisH3.3 and EFTS were the least stable genes in different tissues of H. serrata.

Optimal Number of RGs for Normalization
Though a single and stable RG is su cient for quantifying gene expression, the use of more than one RG for effective normalization of gene expression data is suggested [29]. Based on the geNorm software, the optimal number of RGs needed for normalization was determined by pairwise variation (Vn/n+1). In our data, the all pairwise variation of V2/3 values were lower than 0.15 (Figure 3), which suggested that the combination of the two most stable RGs was optimum for normalization.
Together with RefFinder analysis, GAPDHB and HisH2A were the best combination for normalization in different tissues of H. serrata.

RG Validation
To demonstrate the utility of identi ed stable RGs, four HupA-biosynthetic related genes LDC, MET CYP and CAO were selected in H. serrata. For the purpose of comparison, expression values of target genes were normalized with respect to the most stable gene pair (GAPDHB and HisH2A) and the least stable gene pair (EFTS and HisH3.3) in H. serrata different tissues. When normalized using the most stable genes, the transcription levels of LDC, MET, CYP and CAO were (over 2 -fold) in the tissues of leaf, stem and root were compared, and the expression trend was consistent with that of transcriptome sequencing dates (Figure 4). By contrast, when normalized using the least stable genes, the transcription level of MET and CYP were not up-regulated (less than 2 -fold) in stem and leaf tissues. The transcription level of LDC was down-regulated (0.77 -fold) in stem tissues, and the CAO was down-regulated (0.67 -fold) in leaf tissues. This expression trend was not consistent with that of transcriptome sequencing dates. In all, the expression of the most stable gene pair was more reliable than the least stable gene pair.

Discussion
H. serrata has received extensive concern due to produce biologically active lycopodium alkaloids, especially HupA [7]. HupA was found to possess potent acetylcholine esterase inhibition (AChEI) and had been clinically exploited for the treatment of Alzheimer's disease. More studies are focused on the isolation and identi cation of compounds and endophytic bacteria [31], but little on the transcriptional level. Especially for the selection of optimum reference genes, little research has been reported [35]. So, in present study, we screened and selected the optimal RGs based on full-length transcriptome sequencing and previous researches in different tissues of H. serrata. By the analysis of four softwares and the veri cation of four HupA-biosynthetic genes (LDC, MET, CYP and CAO), we obtained two optimal reference genes GAPDHB and HisH2A for studying HupA-biosynthetic related genes. This study provides suitable normalization for analyzing the expression of HupA-biosynthetic gene. In addition, we found the expression trend of HupA-biosynthetic genes were similar with the trend of HupA content in different tissues of H. serrata. This result will provide the information for further studying the biosynthesis and transportation of HupA. In general, the expression level of RGs should be constant stable in any physical conditions. However, there is no RG suitable to all biological systems. We had to screen the most suitable RGs for studying the HupA-biosynthetic. Based on the CPM value of transcriptome sequencing and reported literatures [32,8,14,5], ten traditional RGs and three new RGs candidate were chosen (Table 1). Especially for three new RGs candidate, they had stable expression in full-length transcriptome sequencing but they were not used as RGs before.
The primer speci city is the primary condition of RT-qPCR. The ideal primers which should cross intron regions to availably avoid genomic contamination in cDNA samples and cannot be set in conservative domain. Firstly, each primer pair was designed except the conserved domains ( Figure S2). Subsequently, the products of each primer pair were detected by agarose gel electrophoresis ( Figure S3) and melting curves ( Figure S4). The results indicated that there were no primer dimers and non-speci c ampli cation for each primer pair. Furthermore, the E value of PCR varied from 90.4-103.6%, and all of the R 2 were greater than 0.990 (Table 2), which were similar to previous literatures [16]. In conclusion, these results indicated each primer pair was speci city and the RT-qPCR assays were highly e cient.
Based on the analysis of four softwares (geNorm, NormFinder, BestKeeper and RefFinder), GAPDHB and HisH2A stood out among the thirteen candidate RGs, and became the best combination for normalization in different tissues of H. serrata (Table 5). Many studies have shown that GAPDH was most often as relatively stable internal RGs in different tissues and under a variety of experimental conditions [6, 24,14]. Although histone and elongation factor were be reported as most stable RGs in other species [24], HisH3.3 and EFTs were the most instable RGs in different tissues of H. serrata. These results suggested that the traditional RG may not suitable for all samples. Interestingly, the performance between the homologous genes was obviously different. HisH3.3 and EFTs were more instable than HisH2A and EF1dt. Similar ndings can be found that the expression level of Actin2/7 was more stable than Actin11 in diverse tissues of soybean [13]. These results also stated clearly that the expression level and stability of RGs from the same gene family may be different in the same samples. Taken together, the results further proved the necessity for screening suitable RGs in different tissues of H. serrata.
A proposed biosynthesis pathway for HupA and related lycopodium alkaloids was reported [18]. However, only two enzymes, LDC and CAO have been proved to participate in the biosynthesis of HupA [33, 4, 28]. Three enzymes RS-norcoclaurine 6-O-methyltransferase (MET) and cytochrome P45072A1 (CYP) [17,33], type III polyketide synthase (PKS) [30], have been described to be possible involvement of the biosynthesis of HupA. In order to verify the accuracy of the stable RGs identi ed in this paper, four HupA-biosynthetic genes, LDC, MET, CYP and CAO were selected. The results showed that when the combination of stable RGs (GAPDHB and HisH2A) was used, the consistent expressions trend of LDC, MET, CYP and CAO were obtained between transcriptome sequencing and RT-qPCR ( Fig. 5and Table 1). Conversely, the use of the most unstable RGs (HisH3.3 and EFTS) may lead to declinational results ( Fig. 5and Table 1). The results further veri ed that GAPDHB and HisH2A were suitable for gene expression normalization, especially for HupA-biosynthetic genes. In addition, we tested the content of HupA in different tissues. The results indicated the content of HupA in root was obviously lower than that in stem or leaf (Fig. 1), which was consistent with the previous study [31,18]. Therefore, in order to protect the wild resources, we suggested picking the aboveground parts instead of uprooting the whole plan when digging H. serrata. In addition, we found the expression trend of HupAbiosynthetic genes were similar with the trend of HupA content in different tissue of H. serrata, which indicated that the biosynthesis of HupA may be in stem and leaf. This result will provide the information for further studying the biosynthesis and transportation of HupA.

Conclusions
In the present study, based on full-length transcriptome sequencing data and the analysis of four softwares, we obtained two optimal reference genes GAPDHB and HisH2A from thirteen candidate reference genes in different tissue of H. serrata. The expression patterns of four HupA-biosynthetic related genes LDC, MET, CYP and CAO further veri ed that GAPDHB and HisH2A were suitable for gene expression normalization. This work provides suitable RGs for the subsequent research of HupA-biosynthetic and transportation in H. serrata. Declarations 32. Wu, Y., Tian, Q., Huang, W., Liu, J., Xia, X., Yang, X. and Mou, H. Identi cation and evaluation of reference genes for quantitative realtime PCR analysis in Passi ora edulis under stem rot condition. Molecular Biology Reports 2020, 47:2951-2962.     The content of Hup A in different tissues of H. serrata. The mean and standard deviation were calculated using the data from three independent biological replicates. Pairwise variation (Vn/n+1) of ten candidate reference genes calculated by geNorm. The cut-off value to determine the optimal number of RGs for qRT-PCR normalization is 0.15.