PCR amplification of nine candidate reference genes.
Nine common reference genes in D. celebensis, alpha-tubulin (Atb), β-actin (Act), 18S ribosomal RNA (18S), glyceraldehyde-3-phosphate dehydrogenase (GAPDH), elongation factor 1-beta (EF-1b), ubiquitin conjugating enzyme (UBC), histone H2A (H2A), TATA-box binding protein (TBP), and succinate dehydrogenase (SDH) were selected as candidate references for the analysis of gene expression under stress conditions made by three representative chemicals (e.g., BPA, B[a]P, Hg), and different ages (24-h, 4-day, 7-day, and 10-day). The information of all nine candidates, such as their gene IDs, gene, size of the PCR products, Tm values, and the median Ct values, is shown in Table 1. In general, the sizes of the qPCR products of all the tested genes ranged from 81bp to 108 bp. The applicability of the primers designed for the qPCR amplification of these nine candidate genes was verified by RT-PCR and sequencing of the amplified products (Supplementary Information, Text S1). The unimodal melting curves indicate that the primers used in this study have high specificity (Text S1). Furthermore, the efficiency values ranged from 90.8–106.6% (Table 1). Therefore, our qRT-PCR results were confirmed to be valid and reliable.
Table 1
Summary of nine reference genes
Gene symbol
|
Gene description
|
GenBank accession no.
|
Primer sequences (5'-3')
|
Amplicon size (bp)
|
Tm (°C)
|
Efficiency (%)
|
Atb
|
Alpha-tubulin
|
MH636293
|
GCATGATTTCCAACACGAC
TACCAGTGAACGAAGGCT
|
97
|
55.2
|
93.92
|
53.9
|
Act
|
beta-actin
|
MH636289
|
CAAGATTGTCGCTCCTCCTG
CATCTGCTGGAAGGTGG
|
84
|
60.5
|
90.79
|
59.5
|
18s
|
18s ribosomal RNA
|
AF144210.1
|
TGGAAGGATTGACAGATTGA
CTTAGTTGGTGGAGCGATTT
|
81
|
54.3
|
101.89
|
56.4
|
GAPDH
|
Glyceraldehyde-3-phosphate dehydrogenase
|
MH636290
|
AACTGTCGCCGCTGTTGA
ATGGAATGTTCTTGGGGTCG
|
86
|
56.1
|
96.31
|
58.4
|
EF-1b
|
Elongation factor 1-beta
|
MH636295
|
CGGCTGTGTCGTTGAAGA
GGCAATGTCCACACTCTG
|
94
|
56.1
|
106.67
|
56.1
|
UBC
|
Ubiqutin Conjugating enzyme
|
MH636294
|
CCTTTTGACGGACCCATAT
TAGTCCAACAGCGAGCAATA
|
104
|
55.2
|
103.89
|
56.4
|
H2A
|
Histone H2A
|
MH636296
|
GAATACCTGGCTGCTGAAG
CAATTGGAGATGACGGGG
|
90
|
56.1
|
101.35
|
57.3
|
TBP
|
TATA-box binding protein
|
MH636292
|
TTTCCTGGCTTAATTTACCGTA
ATCTCTTGTCTGACTTTGGC
|
101
|
54.3
|
96.37
|
56.4
|
SDH
|
Succinate dehydrogenase
|
MH636291
|
CATCGAGTCTCAACAGAAGA
GAGCTCAACCTTTCCAGTT
|
91
|
56.4
|
100.16
|
55.2
|
Expression level of nine reference genes.
As shown in the results and box and whisker plots (Tables 2, 3, and Fig. S1), the Ct value distribution of nine reference genes under different conditions varied from 8.7 to 28.3. The 18S gene encoded the most abundant transcripts in D. celebensis, reaching the threshold fluorescence peak after approximately 10 cycles with the lowest median Ct value (Ct < 10, Table 2, Table 3, and Fig.S1) compared to those of other candidate reference genes during qPCR amplification when the same amounts of total RNA (500ng) were used as templates for reverse transcription reactions. The least abundant transcripts were Atb, which had a Ct value of 21 or higher (21.93 ± 0.38), for chemical treatment. In a test using different ages, UBC had the highest Ct mean values (25.80 ± 0.54). However, with regard to variations in Ct values, Atb genes showed significant differences between all ages of development (ANOVA, p < 0.05).
Table 2
Cycle threshold (Ct) values of nine putative reference genes in 12 mRNA at different conditions (three chemicals and four concentrations) from a monoculture of Diaphanosoma celebensis exposed to chemicals. (Min., minimum; Max., maximum; SD, standard deviation; SE, standard error; CV, coefficient of variation.)
Gene
|
Ct value statistics
|
min
|
max
|
max-min
|
mean
|
SD
|
SE
|
CV
|
Atb
|
21.18
|
22.48
|
1.30
|
21.93
|
0.38
|
0.06
|
0.02
|
Act
|
16.06
|
16.52
|
0.46
|
16.26
|
0.09
|
0.01
|
0.01
|
18S
|
8.70
|
9.61
|
0.91
|
9.19
|
0.21
|
0.03
|
0.02
|
GAPDH
|
17.00
|
18.27
|
1.27
|
17.32
|
0.23
|
0.03
|
0.01
|
EF-1b
|
17.34
|
18.39
|
1.05
|
17.80
|
0.20
|
0.03
|
0.01
|
UBC
|
21.40
|
22.23
|
0.83
|
21.82
|
0.20
|
0.03
|
0.01
|
H2A
|
18.26
|
18.98
|
0.72
|
18.49
|
0.17
|
0.02
|
0.01
|
TBP
|
20.85
|
21.65
|
0.8
|
21.21
|
0.20
|
0.03
|
0.01
|
SDH
|
18.60
|
19.54
|
0.94
|
19.00
|
0.19
|
0.03
|
0.01
|
Table 3
Cycle threshold (Ct) values of nine putative reference genes in 12 mRNA samples from a monoculture of Diaphanosoma celebensis at different ages (Min., minimum; Max., maximum; SD, standard deviation; SE, standard error; CV, coefficient of variation.)
Gene
|
Ct value statistics
|
min
|
max
|
max-min
|
mean
|
SD
|
SE
|
CV
|
Atb
|
18.85
|
28.27
|
9.42
|
23.82
|
3.41
|
0.57
|
0.14
|
Act
|
18.97
|
21.29
|
2.32
|
19.75
|
0.57
|
0.10
|
0.03
|
18S
|
9.59
|
13.43
|
3.84
|
11.18
|
1.14
|
0.19
|
0.10
|
GAPDH
|
19.41
|
21.79
|
2.38
|
20.31
|
0.63
|
0.11
|
0.03
|
EF-1b
|
21.51
|
24.23
|
2.72
|
22.20
|
0.69
|
0.12
|
0.03
|
UBC
|
24.99
|
26.80
|
1.81
|
25.80
|
0.54
|
0.09
|
0.02
|
H2A
|
22.22
|
23.86
|
1.64
|
22.88
|
0.45
|
0.08
|
0.02
|
TBP
|
24.21
|
26.85
|
2.64
|
25.39
|
0.90
|
0.15
|
0.04
|
SDH
|
22.00
|
24.12
|
2.12
|
22.90
|
0.46
|
0.08
|
0.02
|
The stability of candidate reference genes under chemical exposures.
The stability of the reference genes showed different results according to the chemicals (Fig. 1, Fig. S2-S5). Based on the M values evaluated by geNorm, all genes and chemicals showed M values of less than 1, an acceptable value considering suitable reference genes for RT-qPCR normalization in all chemical exposure conditions (Fig. 1, Fig. S2). Act and GAPDH were the best reference genes for samples exposed to B[a]P (Fig. S2A) and BPA (Fig. S2B). Upon exposure to Hg, UBC and TBP were the most suitable reference genes (Fig. S2C). When considering all samples exposed to three different chemicals, the M values of all genes were stable in the following order: H2A, EF-1b, UBC, TBP, Act, 18S, SDH, GAPDH, and Atb (Fig. 1).
NormFinder also showed different results depending on the sample group exposed to the chemical (Fig. S3). Act, Atb, and TBP were the most stable genes for exposure to B[a]P, BPA, and Hg, respectively (Fig. S3A-S3C). The results calculated by NormFinder were somewhat similar to those from geNorm in some representative genes in stability order for each chemical and whole samples (Fig. 1B).
For all cases, including each chemical and all samples, Act was the most stable reference gene in BestKeeper (Fig. 1C and Fig. S4). This method showed similar ranking results to other methods, but the order was observed to be different for Hg exposure (Fig S4C). Upon exposure to Hg, other methods (geNorm and NormFinder) showed that TBP and UBC were suitable, but only the BestKeeper showed Act and H2A as suitable candidates. This result is also the same in the case where all samples were gathered.
RefFinder, a web-based analysis tool, integrates three major computational programs, geNorm, NormFinder, and BestKeeper, and the comparative ΔCt method to comprehensively rank the tested candidate genes. This comprehensive method summarizes H2A and Act as potential reference genes for the chemical treatment of this species (Fig. 1D and Fig. S5).
Interestingly, in BPA exposure, Atb, which showed the most unstable results for the entire sample, was calculated as a stable reference gene as the first (NormFinder) or top ranking in several methods (third stable gene for BestKeeper, second stable gene for RefFinder) only for BPA exposure (Fig. S2-S5). When all three chemicals are considered together, either at the commonly used thresholds (1.5 for M-value for geNorm, SD lower than 1 for BestKeeper) or stability order (RefFinder), H2A and Act have similarly high stability.
The stability of candidate reference genes under biotic conditions: Different ages
Each age group (24 h, 4 d, 7 d and 10 d post hatching) showed different stabilities for the nine reference genes of the brackish water flea D. celebensis (Fig. 2, Fig. S6-9). Interestingly, individuals at post-24h showed TBP to be the most unstable, whereas this gene was relatively stable at other ages of growth (days 4, 7, and 10). Even on day 4, the TBP gene was found to be the most stable in the NormFinder method (Fig. S7). All results showed that Atb was the most unstable, except at 24 h (Fig. 2 and Fig. S6-S9). In particular, the stability of Atb gene expression on day 4 was significantly different from that of the other genes. For example, in geNorm, only Atb had an M value greater than 1, and the next unstable gene, 18S, was 0.5.
NormFinder and BestKeeper also showed that Atb was the most unstable gene, and also showed a large difference in stability value with the next unstable genes, SDH and 18S, respectively (Fig. 2 and Fig S6-S8). NormFinder showed different results from geNorm in the order of stability order: GAPDH, EF-1b, Act, 18S, H2A, SDH, UBC, TBP, and Atb, while the three most stable genes were the same. BestKeeper showed that SDH is the most stable gene when considering all age samples. The stable results of 18S shown in geNorm and NormFinder do not apply to the seven-day individuals as coefficients in BestKeeper of 18S over 1, the threshold of this algorithm. Finally, the integrated approach RefFinder indicates that Act is the most stable gene for covering all ages of D. celebensis (Fig. S8).
Impact of reference genes on real-time qRT-PCR data analysis of target genes.
To validate the selected reference genes using algorithms, the expression of target genes was normalized and compared at different conditions. First, the pairwise variation (V) suggested that two reference genes are required for the normalization of gene expression levels (Fig. 3). In fact, all V values evaluated by geNorm were less than 0.15, which is the criterion for selecting a suitable reference gene number for the normalization of gene expression in chemical treatments. However, a test conducted using individuals of different ages showed that increasing the number of reference genes did not always increase reliability (Fig. 3B). For differential ages, Act and GAPDH were selected using RefFinder as stable reference genes.
The expression levels of GSTmu were compared between the most stable and unstable genes; however, no difference was found in the levels and patterns regardless of stable (H2A) or not (Atb) (Fig. 4). However, studies performed using D. celebenesis of different ages showed significantly different expression patterns of GSTsigma and EcRA depending on which reference gene was used for normalization (Fig. 5 and Fig. S10). In particular, the result using Atb as a reference gene showed superficially different expression patterns of target genes (highly expressed GSTsigma and EcRA at day 4) unlike other cases using stable reference genes (GAPDH, Act, and both genes).