Identification and evaluation of the novel genes for transcript normalization during female gametophyte development in sugarcane

Plant material

The sugarcane (Saccharum spontaneum L.) cultivar Yuetang 91-976 provided by State Key Laboratory for Conservation and Utilization of Subtropical Agro-Bioresources (Guangxi, China) was used for the experiment. When the plants reached the florescence stage, five different stages of the sugarcane female gametophyte (i.e., Archesporial cell (AC), Megaspore mother cell (MMC), Meiosis, Mitosis and Mature) were collected (Fig. S3). The samples were harvested in three biological replicates and quickly frozen with liquid nitrogen before storing at −80 °C.

RNA extraction and RNA-seq

Total RNA was isolated by the Omega Total RNA kit II (R6934–02, USA). The evaluation of RNA quality was performed by gel electrophoresis and 2,000 spectrophotometer assessment at 260 nm (NanoDrop, Thermo Fisher Scientific), and Illumina sequencing was referred to by the method of Zhao et al. (2018). Each sample had three independent biological replicates, and one µg of total RNA was used for library preparation. After sequencing, the adapter sequences and low-quality reads were filtered, and the clean reads were aligned to the sugarcane reference genome by STAR v2.5.0 software. The SourceForge Subread package feature Count v1.5.0 was used for the alignment results for gene quantification.

The cDNA was synthesized using the ThermoScript RT-PCR kit (Life Technologies) in a 20 μL volume reaction under the program: 42 °C for 15 min and 85 °C for 15 s. The qRT-PCR was performed according to the SYBR Premix RT reagent kit system (TaKaRa, Dalian, China).

Candidate gene selection and primer design

Based on the average RPKM values, the CV (i.e. coefficient of variation) values were calculated for each gene under the corresponding experimental stage. This resulted in 50 candidates with the low CV values (Table S1). We considered only the 18 candidates with the highest RPKM values to avoid identifying genes with absent as latent references. The name, ID and sequence information of these genes were obtained from the datas uploaded by Ming Laboratory (https://www.life.illinois.edu/ming/downloads/Spontaneum_genome/). This led to identification of 12 novel selected reference genes (ATPase, FAB2, LOG2, MCB1, MOR1, TIV110, VHA, VLN3, ARM, Expressed, GET4, PP2A, and TIC110) and six traditional genes (EF1α, TUB6, ACT7, ACTIN, GAPDH, and UBC28) (Table 1).

qRT-PCR and statistical analysis

We used the Bio-Rad Real-time PCR (Foster, USA) for qRT-PCR with SYBR Premix Ex Taq plus (TaKaRa, China) in 20 μL volume reaction. The qRT-PCR cycle program applied in this experiment was as follows: 96 °C for 2 min, 40 cycles at 95 °C for 5 s, 60 °C for 34 s, 95 °C for 15 s. Two technical replicates and three independent biological replicates were performed for each reaction, and the negative control reactions without template were performed in parallel. According to the standard curve and formula, the correlation coefficient (R2) and amplification efficiency (E) were computed: E = [10 − (1/slope) − 1] × 100%.

In order to identify whether the expression of each candidate gene was stable, three statistical software packages, geNorm (Vandesompele et al., 2002), NormFinder (Andersen, Jensen & Orntoft, 2004), and BestKeeper (Pfaffl, 2001), were adopted. For geNorm and NormFinder algorithms, the relative expression values of all candidate genes were calculated from the raw Ct values using the formula: E_∆Ct, where ∆Ct = each corresponding Ct value – the minimum Ct value. Then, the resulting values were imported into NormFinder and geNorm to determine the gene expression stability. For BestKeeper analysis, the coefficient variation (CV) and standard deviation (SD) were computed using the raw Ct values. Finally, the comprehensive rankings of the optimal reference genes were acquired by calculating the geometric mean.

Candidate genes selection

Transcriptome analysis across five different developmental stages of sugarcane ovule at AC, MMC, Meiosis, Mitosis and Mature (Fig. S3) resulted in 50 candidate reference genes with the high average read per kilobase per million mapped reads (RPKM) values. SD values across all stages were calculated and then the ratio of SD and Mean Value (i.e., coefficient of variation, CV) were calculated as a reference index. Since lower CV values indicate less deviation, genes with CV values less than 0.6 were picked out from 50 candidate genes. A total of eighteen genes, including twelve novel genes (ATPase, FAB2, LOG2, MCB1, MOR1, TIV110, VHA, VLN3, ARM, Expressed, GET4, PP2A, and TIC110) and six traditional genes (EF1α, TUB6, ACT7, ACTIN, GAPDH, and UBC28) were screened out (Table 1). The CV values of 12 novel genes ranged from 0.06 to 0.09, which were much lower than eight traditional reference genes (0.25–0.56) (Fig. 1A). We further calculated the relative expression level of all candidates across all five experimental stages to obtain more accurate data. Except for EF1α and TUB6, the novel genes showed very stable expression. In contrast, all the traditional housekeeping gene expression levels varied considerably, suggesting that the expression of novel genes may be more stable than conventional housekeeping genes during sugarcane ovule development (Fig. 1). Taken together, the results of RNA-seq showed that most of the new genes are better than traditional housekeeping genes in sugarcane ovule normalization studies.

PCR specificity and amplification efficiency of candidate reference genes

To investigate the transcriptional behaviors of candidate reference genes, RT-PCR and qRT-PCR experiments were adopted to confirm the primer specificity of all candidate genes. In RT-PCR results, the nonspecific amplification band and primer dimer were not observed and no signal was detected in the negative controls of all tested primers (Fig. S1). In the qRT-PCR experiment, all genes had a relatively single melting curve, indicating that the primers used in the experiment had high specificity (Fig. S2). Besides, the standard curve was established with a series of 10-fold diluted cDNA samples and then the correlation coefficient (R²) and amplification efficiency (E) of each gene was calculated. The results showed that E values of all genes varied from 0.97 (ATPase) to 1.04 (VLN3 and TUB6), and the R² of all standard curves ranged from 0.995 (TUB6 and GET4) to 0.999 (GET4, ARM, TIC110 and EF1α), showing their high correlation with their average value (Table 1).

Expression profiles of candidate reference genes

The top eighteen genes were selected as candidate reference genes for qRT-PCR, and the Ct values were examined across all ovule developmental sugarcane samples to detect the transcript abundance and expression stability. The lower the Ct value means the higher the expression level, and vice versa. As showed in Fig. 2, the Ct values of the 10 novel genes ranged between 4.55 and 30.71, while the commonly used housekeeping genes, Actin and ACT7, with the broadest range 6.37 ~ 31.09 (Actin) and 28.36 ~ 7.05 (ACT7), showed the worst expression stability. These results were also consistent with their CV values, 0.54 (Actin) and 0.50 (ACT7), the two highest CV values of all candidates (Table 1). In contrast, the three novel genes FAB2, with the Ct values varied from 26.15 to 30.71, followed by VLN3 (4.97 ~ 23.23) and MOR1 (7.68 ~ 22.71), which showed minimal variation. Therefore, these genes are speculated as more suitable for normalizing sugarcane ovule development.

Stability evaluation by geNorm, BestKeeper and NormFinder analysis

To obtain an accurate result, three standard programs BestKeeper, geNorm and NormFinder were combined to detect the expression stability of eighteen candidate reference genes. For the geNorm algorithm, the average variation of tested genes is defined as M-values, with lower M value and more stability. As shown in Table 2, MCB1 and LOG2 showed the most stable expression level owing to their lowest M values (0.15). While the M value of TUB6 was the highest, indicating that TUB6 is the most variable among these candidates. The NormFinder algorithm is based on a mathematical model that can identify the most suitable reference gene by calculating both the overall variation of all candidates and variation between the given samples (Andersen, Jensen & Orntoft, 2004). The results demonstrated that MOR1, with the stability value of 0.09, followed by FAB2 (0.10), TIC110 (0.10), were the most stable genes. MCB1, which ranked in the first position by geNorm, was ranked in the fourth position by NormFinder. Not surprisingly, TUB6 was also evaluated as the most unreliable gene by NormFinder (Table 2). BestKeeper is also an EXCEL-based tool determining the most stable reference gene by calculating SD and CV. Both SD value and CV value were negatively correlated with gene expression stability, and thus candidate genes with lower SD value and lower CV value exhibit higher normalization reliability (Pfaffl, 2001). As shown in Table 2, FAB2 with the lowest value (CV (%) ± SD = 5.48 ± 1.52) was considered as the most reliable normalization factor, followed by VLN3 (6.28 ± 1.63), TIV110 (11.29 ± 2.68), while TUB6 (20.38 ± 4.32) followed by EF1α (21.69 ± 5.35) were identified as the most unstable normalization factors.

Consensus ranking of candidate reference genes

Owing to different statistical approaches among the three algorithms, the best candidate genes generated by each algorithm were not for their same ranking. To obtain a comprehensive ranking, the geometric mean method was introduced to create consensus results described by Niu et al. (2015). As shown in Table 2, FAB2, MOR1 and MCB1 were considered as the three best genes by comprehensive analysis for sugarcane female gametophyte development.

Moreover, the pairwise variation (V) values were also calculated by geNorm to confirm the minimum number of reference genes required for accurate normalization. According to Vandesompele’s research, when the V values between two adjacent normalization factors (Vn/Vn+1) were higher than the recommended value 0.15, one more reference gene than n is required (Vandesompele et al., 2002). As shown in Fig. 3, the first pairwise variation value is lower than the suggested value 0.15. Therefore, two top-ranked genes are enough for the reliable normalization. The third gene adopted will make no significant difference for accurate normalization across these sugarcane female gametophyte developmental stages. The results of consensus ranking showed FAB2 and MOR1 are the best choices for normalizing transcription during female gametophyte development in sugarcane.

Reference gene validation

To validate the effectiveness of reference genes selected in this study, the expression of two transcripts, CDK and REC8, were analyzed by qRT-PCR analysis across the different ovule stages in sugarcane. CDK encodes a cyclin-dependent protein kinase, playing key roles in the progression of the cell proliferation or cell cycle in eukaryotes (Endo et al., 2012). REC8 is reported to encode a meiotic cohesion protein that could bind sister chromatids and anchor chromosomes to the axis, thus is critical for chromosome segregation during meiosis and eukaryotes gametophyte production (Hsieh et al., 2020). As depicted in Fig. 4, when stable genes FAB2 and/or MOR1 used as the normalization factors, the expression level of REC8 gene increased from AC stage and peaked at the MMC stage, and was reduced at the meiosis stage, followed by an increase at mitosis stage. The trend of relative expression of REC8 and CDK was consistent with RNA-seq data. On the contrary, when the unstable reference genes ATPase and/or TUB6 were selected for validation, the expression profiles of CDK and REC8 showed significant fluctuations, about 3-fold high than that of most genes at MMC stage, and the trend of expression profiles was also changed after the meiosis stage (Fig. 4).