Comprehensive transcriptome analysis of reference genes for fruit development of Euscaphis konishii

Plant material

The fruits of E. konishii have three developmental stages, including the green stage, turning stage, and red stage. These were collected from Qingliu County, Fujian Province, China, from June to September. The samples were provided by Sanming yisheng agricultural and forestry co., LTD. All pericarp samples were separated from the fruit, wrapped in tin foil and then frozen in liquid nitrogen and immediately stored at −80 °C until they were qualified for further analysis. Three biological replicates for each sample were used for RNA extraction.

Establishment of the RNA-Seq database

We sequenced the transcriptome of the E. konishi i pericarp in three developmental stages. The library produced 67.78 G of clean data, and the average clean data in the nine samples were 6.78 G. In this study, high-quality libraries with mapping rate were higher than 79.74% and Q30 values highter than 92.70%. In total, 86,120 unigenes with a mean length of 893.34 bp and N50 length of 1,642 nt were assembled by using Trinity (Grabherr et al., 2011) with min_kmer_cov set to 2 by default and all other parameters set default. The length distribution of all unigenes were shown in Fig. 1. The raw data has been submitted to NCBI, Sequence Read Archive (SRA) submission: SRR9267648 to SRR9267656. We used BLAST software to compare unigene sequences to NR (https://blast.ncbi.nlm.nih.gov/Blast.cgi), Pfam (http://pfam.xfam.org/), GO (http://www.geneontology.org/), KEGG (http://genome.jp/kegg/), Pwiss-Prot (http://uniprot.org/), KOG (http://www.ncbi.nlm.nih.gov/KOG), eggNog (http://eggongdb.embl.de/), and COG (http://www.ncbi.nlm.nih.gov/COG). Transcriptome sequences were used as references to analyse the expression profiles of each sample. The sequencing results of each sample were compared to the reference sequence, and the expression amount of each unigene in different samples was obtained.

RNA-Seq database analysis and primer design

Genes with expression levels lower than 5 FPKM in tanscriptomes would make poor qPCR references due to the difficulties in detecting and quantifying their expression. (Kimmy, Patrick & Joshua, 2017). After their removal, a total of 75,048 genes in E. konishii were evaluated. We calculated the mean expression value, standard deviation, and coefficients of variation based on the raw RNA-seq data. Base on the requirements CV ≤ 0.2 and FC ≤ 0.2, 1,131 genes (Table S1) were obtained by removing overabundant genes with low expression levels. Finally, we selected the eight candidate genes according to genome annotation of each, which was assigned base on the best mach of the alignments using Blast to NR, TrEMBL and KEGG databases. The selecting steps of candidate reference genes were shown in Fig. S1. In order to ensure the accuracy of the reference gene predictions, the coding sequences of eight selected genes was used as queries for BLAST orderly through the TAIR database (http://www.arabidopsis.org/). The sequences with the highest homology with Arabidopsis were shown in Table 1.

Here we obtained the nucleotide sequences of reference genes from our laboratory. The nucleotide sequences of candidatate reference genes were appended in Table S2. According to the nucleotide sequence of candidate reference genes and the design principle of qRT-PCR, forward and reverse primers of all candidate reference genes were designed using Primer Premier 6.0 software and were synthesized by Jinweizhi Biotechnology Co., LTD (Suzhou, China) and purified by PAGE. The primer sequences, amplicon size, anneal temperature were shown in Table 2.

RNA Extraction and cDNA Synthesis

Total RNA was extracted and purified by using an RNAprep Pure Plant Kit (Polysaccharides and Polyphenolics-Rich, Tiangen, Beijing, China), according to the manufacturer’s instructions. Two-hundred nanograms of total RNA of each sample was used as the template. In addition, the cDNA synthesis strand was performed by using a Revert Aid First Strand cDNA Synthesis Kit (Thermo Fisher, Foster City, CA, USA), according to the manufacturer’s instructions and was stored at −80 °C for further experiment.

Candidate reference genes for RT-PCR amplification

The RT-PCR mixture contained 25 µL L of 2 ×EasyTaq^®PCR SuperMix, 0.4 µL L of forward primer (10  µM), 0.4 µL L of reverse primer (10  µM), and 5 µL L of diluted cDNA in a final volume of 50  µL L. PCR conditions: 40 cycles of 3 min at 94 ° C, 94 °C at 30 s, 56 °C at 10 s, 72  °C at 30 s, and 72 °C at 7 min. RNA quality was determined by 1.2% agarose gel electrophoresis.

Candidate reference genes for qRT-PCR amplification

The qRT-PCR mixture contained 3 µL L of diluted cDNA, 5 µL of 2 × SYBR Green Master, 0.4 µL L of forward primer (10  µM), and 0.4 µL L of reverse primer (10 µM), with ddH₂O added to achieve a final volume of 10  µL L. All PCRs were performed using the QuantStudioTM Real-Time PCR System under the following conditions: 40 cycles of 2 min at 95 °C, 5 s at 95 °C and 30 s at 60 °C. The procedure ended with a melt-curve ramping from 60–95 °C. The melting curve was analysed to determine primer specificity. For further confirmation, PCR products were cloned into pGXT vectors respectively and then sequenced by Biosune Biotech (Fuzhou, China).

Statistical analysis

The instrument calculated the Cq value of each sample based on the qRT-PCR experiment. We can analyse the stability of candidate reference genes using GeNorm, NormFinder, and BestKeeper software. Finally, ReFinder was used to calculate the comprehensive ranking of the stability of candidate reference genes based on analysis results.

Validation of the candidate reference genes

Anthocyanin is the key factor affecting fruit colouration as an important plant pigment. In general, fruit colouration is closely related to the content and proportion of anthocyanin. To verify the results of our experiments, the expressions of nine genes related to the anthocyanin synthesis pathway were calculated with the reference genes selected (Fig. S5). Moreover, the EkSS3 gene related to glycometabolism pathway was also used to validate the reliability of reference genes. The details of primer were shown in Table 3. The qRT-PCR experimental method was the same as described above, and the relative expression level was calculated by the 2^−ΔΔct method. Experimental data from three biological replicates were analysed using analysis of variance (ANOVA), followed by Student’s t-test (P < 0.05).

Differentially expressed genes between fruit with varied development

Based on the theree comparisons of green vs. turning, green vs. red, and turning vs. red, we identified a total of 4,804 differentially expressed genes (DEGs). Among them, 2,175 DEGs between green vs. turning, 3,935 DEGs between green vs. red, and 936 DEGs between turning vs. red were detected. Among the 4804 DEGs selected to predict functions by gene ontology (GO) annotation and KEGG pathway analysis. The genes related to the glycometabolism, anthocyanin biosynthesis, flavonoid biosynthesis and carotennoid biosynthesis performed different expression in different stages of fruit ripening.

Reference gene selection and primer design

According to selecting process (Fig. S1), eight candidate reference genes were chosen to perform the gene normalization studies. The sequences of these candidate reference gene were used to design the qRT-PCR primers. The details of eight reference genes were lised in Table 1. The gene models of eight candidate reference refered to the homologous genes of Arabidopsis, and the details were shown in Fig. S2.

Candidate reference genes for RT-PCR amplification

PCR results of the 8 candidate internal reference genes are shown in Fig. S4 with a single band, no primer dimer and nonspecific amplification, which could be used for subsequent qRT-PCR analysis. We sequenced the PCR amplification products of eight genes and blast the correctness of the PCR amplified fragments (Table S3).

Cq values of candidate reference genes

The Cq values for all 8 reference genes are shown in Fig. 2. The Cq values varied from 24.65 (EkGAPDH2) to 35.27 (EkACT7). Moreover, as shown in Fig. 2, EkmMDH2, EkMDH and EkACT7 are more variable than other candidate reference genes. The stability of EkmMDH2, EkMDH and EkACT7 expression is poor.

qRT-PCR analysis

As shown in Fig. S3, the melting curve of the 8 candidate reference genes at different developmental stages only had a single main peak, and the amplification curves between the repeated samples had high repeatability.

Expression stability analysis of reference genes

Evaluating the expression stability of candidate reference genes depended on statistical analysis. Four different statistical software programs (GeNorm, NormFinder, BestKeeper and ReFinder) were commonly used to calculate the variability of the candidate genes expression and determine the one which were the most suitable.

Gene expression stability was determined by the M-value in GeNorm. The gene expression stability increases as the M-value decreases (Wu et al., 2017; Vandesompele et al., 2002). It was considered that candidate reference genes can be used as reference genes when the M-value is less than 1.0; the result of GeNorm revealed that EkUBQ1 was the most stable gene with the lowest M-value, followed by EkUBC23, EkCYP38 and EkGAPDH2. However, EkACT7, EkmMDH2, EKMDH and EkTUA3 were unsuitable as reference genes, with a value greater than 1.0. Expression stability values analysed by GeNorm are shown in Fig. 3A.

NormFinder software ranks all reference gene candidates based on intra- and inter-group variation and combines results into a stability value for each candidate reference gene. A better stable gene expression is indicated by a lower stability value (Andersen, Jensen & Ørntoft, 2004). The analysis results of NormFinder shown in Fig. 3B revealed that EkGAPDH2 was the most stable gene for all samples.

BestKeeper software determines the most stably expressed genes based on the standard deviation (SD). The lower the SD, the greater the reference gene expression stability will be (Pfaffl et al., 2004). As the analysis results in Fig. 3C show, the SD values of EkACT7, EkmMDH2, EkMDH, EkTUA3 and EkGAPDH2 are greater than 1 and are considered unacceptable as reference genes, according to the selection criteria of BestKeeper software. EkUBC23, EkUBQ1 and EkCYP38 are stable genes with low SD values.

ReFinder software was used to rank the stability obtained by the analysis of GeNorm, NormFinder and BestKeeper on the comprehensive index. The stability of the reference gene expression is directly related to the size of the index. According to the analysis result of ReFinder shown in Fig. 3D, the stability ranking obtained by ReFinder software is as follows: EkUBC23, EkCYP38, EkGAPDH2, EkUBQ1, EkMDH, EkTUA3, EkmMDH, and EkACT7. The results showed that the stability of EkUBC23, EkCYP38 and EkGAPDH2 was better, and the low expression reference genes (EkUBC23 and EkCYP38) were favourable for quantifying their targets, while the high expression reference gene (EkGAPDH2) was beneficial for quantifying its target.

The suitability of the reference gene

According to the results of the four algorithms (GeNorm, NormFinder, BestKeeper and ReFinder), EkUBC23, EkCYP38 and EkGAPDH2 performed more stably. EkUBC23 and EkCYP38 showed similar expression levels during fruit developmental stages. However, EkUBC23 was used to normalize the expression due to its greater stability. Relative expression levels were normalized using the low expression reference gene (EkUBC23) and high expression reference gene (EkGAPDH2). The results showed that there were some differences between the expression levels of 9 genes related to the anthocyanin synthesis pathway and transcriptome sequencing (Fig. 4). When EkUBC23 was used as the reference gene, the results had no significant difference. When EkGAPDH2 was used for normalizing low expression target genes [c57877.graph_c0( DFR), c59825.graph_c0(CHS) and c69862.graph_c1(UFGT) ], the results were different from the expression of transcriptome. However, EkGAPDH2 has better accuracy than EkUBC23 when it was used for normalizing high expression target genes [c72659.graph_c0(CHS) and c72737.graph_c0(CHI) ]. Therefore, we suggest that the stability and expression of reference genes should be considered as important selection conditions.

To ensure the reliabilty of result, We chosen EkSS3 gene related to glycometabolism pathway as target gene. Relative expression levels were normalized using two most stable reference genes (EkGAPDH2 and EkUBC23). As shown in Fig. S6, when normalized using EkGAPDH2 and EkUBC23 as reference gene respectively, the results had no significant difference. Compared with the green stage, EkSS3 expression was clearly downregulated at turning stage and red stage. The EkGAPDH2 and EkUBC23 performed well when normalizing target gene in different pathway.