Selection and Evaluation of Reference Genes for Reverse Transcription-Quantitative PCR Expression Studies in a Thermophilic Bacterium Grown under Different Culture Conditions

The phylum Deinococcus-Thermus is a deeply-branching lineage of bacteria widely recognized as one of the most extremophilic. Members of the Thermus genus are of major interest due to both their bioremediation and biotechnology potentials. However, the molecular mechanisms associated with these key metabolic pathways remain unknown. Reverse-transcription quantitative PCR (RT-qPCR) is a high-throughput means of studying the expression of a large suite of genes over time and under different conditions. The selection of a stably-expressed reference gene is critical when using relative quantification methods, as target gene expression is normalized to expression of the reference gene. However, little information exists as to reference gene selection in extremophiles. This study evaluated 11 candidate reference genes for use with the thermophile Thermus scotoductus when grown under different culture conditions. Based on the combined stability values from BestKeeper and NormFinder software packages, the following are the most appropriate reference genes when comparing: (1) aerobic and anaerobic growth: TSC_c19900, polA2, gyrA, gyrB; (2) anaerobic growth with varied electron acceptors: TSC_c19900, infA, pfk, gyrA, gyrB; (3) aerobic growth with different heating methods: gyrA, gap, gyrB; (4) all conditions mentioned above: gap, gyrA, gyrB. The commonly-employed rpoC does not serve as a reliable reference gene in thermophiles, due to its expression instability across all culture conditions tested here. As extremophiles exhibit a tendency for polyploidy, absolute quantification was employed to determine the ratio of transcript to gene copy number in a subset of the genes. A strong negative correlation was found to exist between ratio and threshold cycle (CT) values, demonstrating that CT changes reflect transcript copy number, and not gene copy number, fluctuations. Even with the potential for polyploidy in extremophiles, the results obtained via absolute quantification indicate that relative quantification is appropriate for RT-qPCR studies with this thermophile.

Introduction evident that commonly employed reference genes-thought to be constitutively expressed and minimally regulated-in fact display considerable variance across time, cell type, or experimental treatment [16]. This may be attributed in part to the fact that some housekeeping proteins are involved in functions other than basal cell metabolism [17,18] or that housekeeping metabolism is a very dynamic process that is efficient at adapting to different growth conditions [19]. Additionally, the use of multiple reference genes is now widely advocated, as the use of a single reference gene has been shown to result in relatively large errors [15].
Limited studies exist that focus on reference gene selection in extremophiles, with only a single dedicated study performed in recent years for the acidophile Acidithiobacillus ferrooxidans [28]. Due to their unique physiologies, suitable reference genes in extremophiles may differ from those identified in mesophiles. For example, temperature stress in multiple mesophilic species identified elongation factor TUF 6-phosphogluconate dehydrogenase (6PGDH) as most stably expressed [29]-yet the temperature to which the species were exposed (40)(41)(42).5°C) is below the limit at which thermophiles thrive. Thus, the selection and evaluation of appropriate reference genes is a prerequisite when utilizing expression profiling to identify key pathways and mechanisms in Thermus sp. grown under different environmental conditions.
One additional issue that must be addressed for gene expression profiling under different growth conditions or over time in both extremophiles and mesophiles is that of genome copy number. In general, bacteria are thought of as monoploid (i.e. harboring only a single copy of their chromosome). However, the number of polyploidy "exceptions" may now exceed the monoploidy "rule." [30]. A survey of 11 species of proteobacteria demonstrated that only three of the species were monoploid, with the gamma-proteobacterium Pseudomonas putida containing 20 copies of its genome [31]. The extremely thermophilic bacterium Thermus thermophiles HB8 has been shown to display polyploidy, harboring multiple genomic copies in a cell [32]. Similarly, Deinococcus radiodurans, an exceptionally radioresistant bacterium closely related to Thermus sp., has also been shown to be a polyploidy organism, possessing four genome copies per cell in the stationary phase and up to 10 copies per cell during exponential growth [33,34]. Even species typically classified as monoploid may display polyploidy, as genome copy number may vary depending on growth rate: cells that are growing more quickly have multiple initiation events, and therefore multiple genomes per cell. This has been shown to be the case with P. aeruginosa, which contains 2-3 copies of the genome per cell during exponential growth, and 1-2 copies during stationary phase [35]. Genome copy number in the model bacterium Escherichia coli, long thought to be monoploid, has also been shown to increase to 6-7 copies as a function of growth rate, while the cyanobacterium Synechocystis possesses 218 genome copies in the exponential phase and 58 during stationary phase [30].
It is necessary to determine whether changes in gene expression are an accurate portrayal of transcript fluctuations, or whether changes in transcript abundance are derived from gene copy number fluctuation due to potential polyploidy. This can be accomplished via absolute quantification of both transcript and gene copy numbers standardized to a unit value such as cell number, volume of sample, or mass. This transcript-to-gene copy number ratio provides a means to determine whether transcript variation is derived from changes in gene or transcript copy number. We previously used this approach to identify new reference genes and confirm that relative quantification was an appropriate method for assessing gene stability with fluctuations in gene copy number with the fungus Neurospora crassa, which possesses multinucleated cells at various stages during growth [36]. Little to no information exists as to the appropriate reference genes to use in Thermus sp. gene expression studies. Additionally, the issue of gene expression changes as a result of fluxing genome copy number in extremophiles has not been addressed. The objective of the current study was to examine the potential of a suite of genes for use as reference genes over time and/or under different culture conditions in T. scotoductus. The expression stabilities and appropriate rank of these 11 genes were evaluated using the Best-Keeper [37] and NormFinder [38] software tools. Transcript and gene copy numbers were then examined using absolute quantification in a subset of these genes in order to determine whether instances of expression instability were a result of changes in transcript or gene copy number.

Materials and Methods Organism Maintenance
Thermus sp. SA-01 (ATCC 700910) was obtained from the American Type Culture Collection. Lyophilized stocks were initially reconstituted in Castenholtz TYE (ATCC medium 416) per ATCC recommendations. Stock cultures were stored at -80°C in 20% glycerol. Cultures were routinely cultured under aerobic conditions at 65°C by inoculation of a single colony into TYG broth (5 g tryptone, 3 g yeast extract, 1 g glucose per liter per the methods of [7]) followed by 1% inoculation into the desired culture medium. For anaerobic culture experiments, the aerobic TYG culture was passed through one sub-culture consisting of a defined basal medium [6] amended with 30 mM glucose as the carbon source and 10 mM sodium nitrate as the terminal electron acceptor.

Culture Conditions and Cell Harvesting
Thermus sp. was incubated under a range of culture conditions in order to assess growth and test the expression stability of a suite of potential reference genes over time or when in the presence of different electron donors and acceptors. All experiments were conducted at 65°C. Table 1 summarizes the culture conditions. A minimum of three of biological replicates were performed under each condition. Briefly, reference gene expression stability was assessed (1) over time in a specific culture condition or (2) at specific time point(s) between different growth conditions. Anaerobic growth experiments were carried out in 50 mL Hungate tubes (Bellco Glass, Inc.) filled with 10 mL of the appropriate medium (described below). Filter-sterilized carbon sources and terminal electron acceptors were added to the appropriate concentrations, and the tubes sealed with black butyl rubber stoppers under an N 2 /CO 2 atmosphere. Aerobic growth experiments were carried out in 50 mL glass tubes filled with 10 mL of TYG. Both anaerobic and aerobic cultures were grown under mild agitation (70 rpm). Growth was also compared between two different heating methods: incubation in an oven (i.e. traditional incubator) versus heating within a synthetic microwave (CEM) set to a power of 300W. These experiments utilized 60 mL TYG in 100 mL Teflon digestion vessels (CEM) with no agitation.
At each time point, cells were collected for DNA and RNA extraction. For RNA extraction, 4 mL of sample was centrifuged for 3 min at 12,000 X g, the pellet washed in 500 μL RNAProtect Bacteria Reagent (Qiagen), centrifuged for an additional 3 min at 12,000 X g, and re-suspended in a final volume of 40 μL RNAProtect Bacterial Reagent. Samples were stored at -80°C until total RNA extraction.

DNA Extraction
Genomic DNA was extracted from biological triplicates using the DNeasy Blood and Tissue Kit (Qiagen, Valencia, CA) following the manufacturer's instructions for pre-treatment of Gramnegative bacteria. The initial Proteinase K incubation was performed for 2 hours at 56°C, with frequent vortexing. Samples were eluted in a final volume of 100 μL Buffer AE and stored at -20°C until use. DNA concentration and purity were assessed using the Nanodrop-2000 (Thermo Fisher Scientific, Wilmington, DE). DNA was diluted 1:5 in nuclease-free water prior to use.

RNA Isolation and cDNA synthesis
An optimized protocol was developed for total RNA extraction from Thermus sp. based on the general methodology of the RNeasy Mini Kit (Qiagen, Valencia, CA). Briefly, 180 μL of lysozyme (15 mg/mL) re-suspended in TE Buffer (30 mM Tris-HCl, 1 mM EDTA, pH 8) plus 20 μL Proteinase K (Qiagen) was added to each cell pellet (previously re-suspended in 40 μL RNAProtect) and samples incubated on a shaking platform (125 rpm) for 25 min at room temperature. Seven hundred μL Buffer RLT containing β-mercaptoethanol was then added to each sample. The remaining steps followed the manufacturer's protocol. The optional on-column DNase digestion was extended to 1 h. An additional spin at top speed (21,130 X g) for 2 min was included immediately prior to elution in 40 μL RNase-free water. Total RNA concentration and purity was measured using the Nanodrop-2000.
Total RNA was converted to cDNA using the High-Capacity RNA-to-cDNA kit (Applied Biosystems, Foster City, CA) using a previously-optimized protocol [36] with slight modifications. Each reaction contained: 2 μL 20x enzyme mix, 2x RT buffer, a volume equivalent to 1 μg total RNA, and brought to a final volume of 40 μL with nuclease-free water. To check for the presence of DNA, a no-RT control, consisting of all components except the RT enzyme, was performed for each sample. Reactions were incubated in a 2720 thermocycler (Applied Biosystems, Foster City, CA) under the following conditions: 37°C for 60 min followed by heating to 95°C for 5 min. Samples were diluted 1:5 in nuclease-free water prior to use in qPCR.

Candidate reference gene primer design and validation
Quantitative PCR assays were developed for a suite of candidate reference genes (Table 2) based on the SYBR Green chemistry [39]. For each assay, primers were designed that targeted a ca. 100-130 bp region of the coding sequence of each gene. Sequences for each gene were obtained from the annotated Thermus SA-01 genome (NCBI accession number CP001962) available in GenBank. Primers were designed using the Primer3 software (http://bioinfo.ut.ee/ primer3/), with an optimal annealing temperature of 60°C. Primers were examined for potential secondary structures using Mfold [40]. Self and cross dimer formation was examined for each of the primer sets using the Operon oligonucleotide analysis tool (https://www.operon. com/oligos/toolkit.php), while the potential for self-complementarity was examined using an online oligonucleotide properties calculator (http://www.basic.northwestern.edu/biotools/ oligocalc.html). Primers were obtained from Invitrogen (Carlsbad, CA). Primer sets were optimized over concentrations spanning 150-300 nM, with an annealing/extension temperature of 60°C. Product specificity was determined via melt curve analysis. The PCR efficiency of each assay across all culture conditions was determined using LinRegPCR [41].

Construction of external standards for absolute quantification
Absolute quantification was employed to determine the ratio of transcript to gene copies per mL as an additional means of examining reference gene stability. The expression of a subset of the candidate reference genes was also assessed via absolute quantification of the ratio of transcripts to gene copies per mL of sample. This was achieved through the creation of external calibration curves for a subset of the genes. Standards derived from PCR products were constructed for the gyrB, polA2, and rpoC genes, each of which occurs as a single copy within the T. scotoductus genome, and used for absolute quantification of both transcript and gene copy number. Standards were created via amplification from genomic DNA using gene-specific primers (S1 Table). Linearized PCR products were used rather than plasmid-or RNA-based standards as it has been shown that external calibration curves constructed with DNA rather than RNA serve as better models for mRNA quantification due to higher sensitivity, increased quantification range, and higher reproducibility and stability than the RNA calibration curve [42]. Additionally, plasmid-based standards have been shown to result in large overestimation compared to linear DNA in qPCR [43]. PCRs were conducted in 25 μL reaction volumes with the PCR Reagent System (Invitrogen, Carlsbad, CA) in a 2720 thermocycler (Applied Biosystems, Foster City, CA) and contained (as final concentration): 1x PCR buffer (minus Mg), 2 mM MgCl 2 , 0.25 mM each dNTP, 200 nM each forward and reverse primer, 0.5 U Taq polymerase, and brought to a final volume of 25 μL with nuclease-free water. Thermocycling conditions consisted of 94°C, 3 min; 35 cycles of 94°C for 45 s, 55°C for 30 s, 72°C for 45 s, and a final extension at 72°C for 7 min. PCR products were visualized by gel electrophoresis to confirm product amplification of the proper size. The resulting amplicons were purified with the QiaQuick PCR Purification kit (Qiagen), eluted in nuclease-free water, and sequenced with the gene-specific primers by GeneWiz, Inc. (Germantown, MD). Sequence identity was confirmed using the standard nucleotide BLAST function on the NCBI website. The number of copies per μL was calculated using the equation: X g μL -1 DNA/(PCR amplicon x 660) x 6.022 x 10 23 . Serial dilutions for each assay were prepared as described previously [36], so that the final linear range of each assay spanned from 1 x 10 2 to 1 x 10 9 copies. The standard curve for each assay was obtained by plotting the cycle at which fluorescence for that sample crossed the threshold value (cycle threshold, C T ) against the calculated copy number (S1-S3 Figs). The slope, y-intercept, r 2 value, and PCR efficiency of each assay can be found in the Supporting Information (S1 Table).

Quantitative PCR of candidate reference genes
The Fast SYBR Green Master Mix (Applied Biosystems, Foster City, CA) was used for all reference gene qPCR assays. Each reaction contained: 10 μL 2x Fast SYBR Green master mix, 150-300 nM each forward and reverse primer, 5 ng (2 μL) template cDNA, and brought to a final volume of 20 μL with nuclease-free water. Reactions were assembled into 384-well plates using the QIAgility (Qiagen) automated pipetting system. Reactions were performed on the Vii 7 Real Time PCR System with the 384-well block format (Applied Biosystems, Foster City, CA). The following protocol was used for all assays: an initial 20 s incubation at 95°C, followed by 40 cycles of 95°C for 1 s and 60°C for 20 s, followed by a melt curve analysis of 95°C for 15 s, 60°C for 1 min, and 95°C for 15 s to determine product specificity. All qPCR reactions were performed in triplicate using Applied Biosystems MicroAmp Fast 384-well reaction plates sealed with MicroAmp optical adhesive film. No-RT controls were performed with all samples using the rpoC assay (S1 Data). No-template controls were also included in each amplification run to monitor for contamination. Reactions were recorded and analyzed using the Applied Biosystems ViiA 7 System software. The raw threshold cycle (C T ) values and intra-assay variations are provided in S1 Data.

Reference Gene Data Analysis
Three general types of culture conditions were used to determine the appropriateness of a suite of genes for use as reference genes: (1) aerobic growth in a rich (TYG) medium (10 mL culture volume); (2) anaerobic growth in basal medium supplemented with various electron donors and acceptors (10 mL culture volume); and (3) aerobic growth in an incubator versus microwave heating (60 mL culture volume). Samples were collected over a ca. 30-h time span for small-volume cultures, and 48 h for large-volume oven/microwave cultures. The expression levels of the candidate reference genes were initially determined from their C T values. The BestKeeper v. 1 and NormFinder software programs were used to determine the expression stability of all 11 candidate reference genes under each of the different culture conditions and among all conditions. For data input into NormFinder, C T s of each gene at each time point were converted to relative quantities as described previously [36] using the equation Q = E ΔC T , where ΔC T = C T highest abundant sample-C T sample, E = PCR efficiency, and Q = relative quantity. It should be noted that in NormFinder, a low stability value indicates a more stable gene, while in BestKeeper, the coefficient of correlation should be as close to one as possible for a stable gene.
To determine whether instability was due to changes in gene expression or the result of multiple gene copies, indicative of polyploidy, absolute quantification was employed for a subset of the genes: gyrB, polA2, and rpoC. Both the number of transcripts (RNA) and gene copies (DNA) per mL was determined for each gene and used to calculate the transcript: gene copy number ratio under each of the culture conditions. The number of transcripts and gene copies per mL were calculated via external calibration curves specific for each target. Transcript or gene copy numbers from total RNA or DNA were calculated automatically from the regression lines by the ViiA 7 System software (contained in S1 Data). Gene copies per mL were calculated using an equation described previously [44], in which gene copies per mL sample = (gene copies per reaction mix x volume of DNA [μL])/ (DNA per reaction mix [μL] x mL sample used) and adjusted to account for the 1:5 initial dilution factor (volume for each sample listed in S2 Data). Each reaction utilized 2 μL of template from a 100 μL elution. Transcript copies per mL were calculated using an equation described previously [36], with a slight modification to adjust for a starting concentration of 1 μg total RNA and no initial 1:5 dilution of total RNA prior to the RT reaction. Therefore, the following equation was used to calculate transcript copies per mL: [(copies per μL per reaction mix)/(volume μL total RNA added to RT reaction/RT reaction volume μL x cDNA dilution factor] x vol total RNA eluted (μL) (S2 Data). Each reaction utilized 2 μL template from a 40 μL elution. The transcript to gene ratio was determined by dividing transcript copies per mL by gene copies per mL as described previously [36]. Ratio and C T values of all three genes under all conditions were first examined for normal distribution (Shapiro-Wilkes) and equal variance (Brown-Forsythe) (SigmaPlot v.13) under each of the conditions compared (i.e. small-volume aerobic/anaerobic conditions; terminal electron acceptors; and large-volume aerobic oven/microwave incubation). If both ratio and C T data passed these tests, the Pearson Product Moment correlation test was applied to examine the strength of correlation between the two. If one of the data sets failed the test for normality, the Spearman rank order correlation test was applied. T-tests or one way ANOVAs were used to test for significant differences among the ratios of biological triplicates at individual time points between culture conditions (for example, 6 h time point comparing oven versus microwave incubation) or across time under each specific condition (for example, one-way ANOVA to all ratios at all time points of cells grown aerobically in conventional oven). As with correlation analysis, if data did not pass the tests for normality and equal variance, the corresponding nonparametric tests were applied (Mann-Whitney U and Kruskall-Wallace, respectively).

Results and Discussion
General characteristics of candidate reference gene quantitative PCR assays Numerous studies, across different organisms, have shown that the typical reference genes are regulated and vary under different conditions [14,45]. Studies performed thus far have failed to identify a single gene that retains sufficient overall expression stability suitable for all contexts, and that the best candidates differ among environmental conditions [46]. Additionally, while reference gene data exist for the well-studied extremophiles T. thermophilis and D. radiodurans, minimal information exists for other extremophilic species. In the present study, 11 candidate reference genes were selected for evaluation with T. scotoductus, a thermophile shown to be capable of metal reduction and anaerobic growth in addition to aerobic growth in rich medium. Gene selection was based on a combination of parameters including their role within the cell and their use in the literature. Additionally, genes were selected from different areas of the genome (based on annotation of strain SA-01) to minimize the chance of transcriptional coupling affecting the results. For example, rpoC commonly serves as a reference gene in bacterial expression studies (often times without proper validation); a glycosyl hydrolase served as a reference gene for D. radiodurans in transcriptional studies upon exposure to radiation [47]; glyceraldehyde 3-phosphate dehydrogenase (gap) served as a reference gene in this same organism when studying DNA strand break repair [48]; rpoC, the chaperone dnaK, and genes coding for translation initiation factors were used in prior expression studies concerning T. thermophilis bacteriophage [49]; and rpoC and the gyrase subunit A (gyrA) were employed in transcriptome studies with the extremophile Acidithiobacillus ferrooxidans [28]. gyrA has also been evaluated as a reference gene with other bacterial species, with conflicting results (i.e. inappropriate for use with Staphylococcus aureus [50], appropriate with the plant pathogen Pectobacterium atrosepticum [51]).
The compliance of the RT-qPCR assays with the MIQE (Minimum Information for Publication of Quantitative Real-Time PCR Experiments [52]) guidelines is shown in the MIQE checklist (S2 Table). All assays were specific for their intended gene target, as indicated by a single peak in the melt curve analysis (S2 Fig) and possessed similar PCR efficiencies ( Table 2). The average C T values, when assimilating the data from all 11 genes, ranged from ca. 16

Reference Gene Stability during Aerobic and Anaerobic Growth Conditions
The expression stability of candidate reference genes was initially evaluated based on calculated variations (SD, standard deviation of the C T , and CV, the coefficient of variance expressed as a percentage on the C T level) using BestKeeper. In general, a lower SD equates to greater stability, and any gene with the SD > 1.00 is considered inconsistent [37]. In addition to stability (SD and CV) values, additional key parameters include the Pearson correlation coefficient (r), derived from pair-wise correlation analyses that estimate inter-gene relationships among all reference genes, and their associated probability (p) values. The comprehensive results of the BestKeeper analysis are provided in S3 Data. The candidate reference genes were then ranked to determine the best combination of genes under each condition and across all conditions using the program NormFinder. One of the valuable functions of this program, in addition to assessing overall expression variation, is that the samples can be assigned to different groups, so that the intra-and intergroup expression variation of the candidate genes can be assessed. This allows for the identification of specific conditions under which each gene displays greatest variance. When comparing reference gene expression across time among aerobic growth in a rich medium (TYG) and anaerobic growth in basal medium supplemented with glucose or lactate as the carbon source (with nitrate as terminal electron acceptor in both), stability values as determined via BestKeeper ranged from 0.80-1.54. The glycosyl hydrolase TSC_c19900 (hereafter referred to as "19900") ranked as the most stable across time and all conditions, followed by polA2 (0.90), gyrA (0.99), gyrB (1.00), and gap (1.02). These were the only genes with a SD 1.00. A gene commonly employed as reference gene, rpoC, was determined to be the least stable among these conditions (SD = 1.54), along with the glycosyl hydrolase TSC_c02940 (hereafter referred to as "02940") (SD = 1.54) ( Table 3). Growth conditions were analyzed three different ways with NormFinder: (1) groups assigned by medium type (glucose/nitrate, lactate/nitrate, or TYG); (2) groups assigned by time (early, mid, and late growth); and (3)  groups assigned by both medium and time ( Table 3). The most appropriate gene when samples were assigned by medium type was polA2 (0.137), with the best combination of genes identified as gap and polA2 (0.144) (S4 Data). With regards to NormFinder analysis, potential reference genes should display an inter-group variance as close to zero as possible [38]. However, it is important to note that the simultaneous use of genes with positive and negative inter-group values would yield lower stability values than sets involving genes with the same sign due to compensating effects [53]. All genes displayed positive and negative inter-group values, with a summation totaling zero (S3 Table). When samples were grouped by time, polA2 was also identified as the most stable (0.059), with the best combination being gap and polA2 (0.047) (S4 Data). All genes also displayed positive and negative inter-group values by this grouping mechanism, with summation of all three variables totaling zero (S3 Table). PolA2 was also identified as the gene with the lowest stability value (0.208) when grouping by both time and medium, with the best combination identified as gap and polA2 (0.195) (S4 Data). Gly was identified as the least stable gene. As with the other groupings, grouping by both time and medium resulted in a summation of zero when totaling the intergroup values for all 9 variables (S3 Table). Analysis with no grouping also identified polA2 as most stable (0.116), followed by gap (0.197). With samples not defined by groups, NormFinder identified gly as the least stable (0.917), followed by 02940. Reference Gene Selection for RT-qPCR with Thermus sp.

Reference Gene Stability during Anaerobic Growth with Varied Electron Acceptors
Under conditions of anaerobic growth with different terminal electron acceptors, the most stable genes as determined via BestKeeper were 19900, infA, and pfk (Table 4). All three displayed very low and very similar SD values of 0.68, 0.70, and 0.72, respectively. As the PCR efficiencies of these assays are also very similar (1.84, 1.83, and 1.81, respectively [ Table 2]), any of these three could be used alone or in combination as appropriate reference gene in the short-term (i.e. 24 h growth) expression analyses under anaerobic growth with different terminal electron acceptors. Following these three genes were the gyrase subunits, gyrA (0.81) and gyrB (0.86). Similar to the results obtained when comparing aerobic and anaerobic growth over time, rpoC would serve as the least appropriate reference gene under these conditions as well, as it was the least stable by a relatively large margin (SD = 1.89) ( Table 4) among the set. NormFinder identified 02940 as having the lowest stability value (0.224), with the best combination being 02940 and 19900 (0.225) (S4 Data). In agreement with the BestKeeper analysis, these were followed by gyrB and infA. All genes displayed a cumulative intergroup variation of zero, as all displayed both a positive and negative inter-group value. When considering the intragroup variation of cells incubated with iron as the terminal electron acceptor, gap, 02940, pfk, and dnaK displayed the lowest values, indicative of stability among samples under this condition, while 19900, infA, gyrB, and polA2 exhibited the greatest stability among samples grown in the controls in which nitrate served as the terminal electron acceptor (S4 Table).

Reference Gene Stability during Oven or Microwave Incubation
When comparing heating methods (oven versus microwave) at 65°C, BestKeeper analysis identified gyrA (0.75), gap (0.80) and gyrB (0.84) as the most stable genes. The SD of multiple genes were equal to or slightly greater than the recommended cut-off of 1.00 (dnaK, infA, polA2) ( Table 5). As with the results obtained under varying anaerobic growth conditions as well as aerobic growth, rpoC was the least stable gene (SD = 1.83), along with 02940 (Table 5). Oven and microwave-grown samples were then analyzed by different groupings in NormFinder: (1) heating method (oven or microwave) (2) time (3) both variables combined and (4) no groupings. The most appropriate reference gene based on heating method was polA2 (0.071), with polA2 and 19900 as the best combination (0.120) ( Table 5 and S4 Data). In accordance with these results, polA2 displayed the lowest intragroup variation among all assays under either heating method (MW = 0.001, oven = 0.002), followed by infA (0.003), dnaK (0.014), and  Table). PolA2 (0.115) was also the most stable gene across time, with infA and polA2 as best combination (0.088) (S4 Data). Supporting these individual grouping analyses, polA2 was also the most stable based on combined heating method and time (0.138), with polA2 and dnaK as best combination (0.141) (S4 Data). The intra-and intergroup values when grouping by heating type; time; and a combination of heating type and time can be found in the Supplementary Information S5 Table. Reference gene stability across all conditions Assimilating the results of both BestKeeper and NormFinder under the various culture conditions, different genes were identified as the most stable under different conditions. While gyrB and 19900 both ranked among the top four in two out of the three treatment comparisons (gyrB, aerobic versus anaerobic growth and oven versus microwave heating; 19900, aerobic versus anaerobic growth, anaerobic growth comparing terminal electron acceptors) there was not a single gene that was consistently identified as the most stable among all conditions. Therefore, in order to accurately identify the most stable genes under all conditions, all data were combined and analyzed with both BestKeeper and NormFinder. The results of BestKeeper revealed that 3 of the 11 genes displayed acceptable stability across all conditions (reflected by a SD <1.00): gap, gyrA, and gyrB; additionally, infA, polA2, and 19900 were just slightly above the cut-off of 1 (all <1.20) (Table 6). Conversely, the gene most frequently used for normalization, rpoC, displayed the greatest variance. The most appropriate gene for use across all conditions as identified via NormFinder was gap, with the best combination consisting of gap and infA, with a combined stability value of 0.206 (Table 6 and S4 Data). PolA2, gyrA, and gyrB were the next most-stable genes among the 11 according to NormFinder. In keeping with the positive and negative values obtained for the intergroup variation and the resulting summation of zero, all genes also displayed a summation of zero when combining all conditions (S6 Table).

Absolute quantification of transcript and gene copy numbers
Other bacterial species closely related to T. scotoductus (T. thermophilis, D. radiodurans) have been shown to exhibit polyploidy. Additionally, growth stage and/or environmental factors have been shown to stimulate polyploidy in bacterial species not typically associated with this trait (i.e. E. coli, P. aeruginosa). Based on this, expression fluctuations in a subset of the potential reference genes were examined with absolute quantification in order to investigate whether fluctuations in C T s were indicative of changes in transcript abundance, or whether C T fluctuations were an indication of changes in gene copy number. This would be the case if cells were transitioning between mono-and polyploidy. The three genes were selected based on a combination of factors, as they play different roles within the cell and ranged in their stability values under the various growth conditions. External, PCR-amplicon-based standards were created and used to calculate both transcript and gene copy numbers per mL of sample in order to determine the ratio of transcript copies to gene copies under the various culture conditions (representative image in Fig 2, plots of transcript and gene copy number for all genes under all conditions S5-S13 Figs). This transcript: gene copy number ratio indicated whether C T fluctuations were due to changes in transcript abundance, indicative of changes in gene expression, or changes in gene copy number, such as would occur during mono/polyploidy transitions. As described in a previous study, the Pearson Product Correlation test was applied to the ratio and corresponding C T of each time point for each gene to test the hypothesis that a negative correlation exists between ratio and C T [36].
With the exception of gyrB across small-volume aerobic/anaerobic conditions ( Table 7, "gyrB Growth"), statistically significant negative correlations were found to exist between the transcript: gene ratio and C T values among the three genes under all conditions (Table 7). These correlations provide strong evidence that C T changes are an accurate depiction of transcript copy number, and not gene copy number, fluctuations. Plotting both the transcript: gene copy number ratio and C T values over time illustrates that in general, C T changes reflect fluxes in transcript copy number (Figs 3-5). For example, when grown aerobically in an oven, the gyrB C T and ratio values showed an inverse trend over time (Fig 3E). If this were due to changes in gene (DNA) copy number, the ratio pattern would mirror that of the C T s. However, an increase in ratio is strongly indicative of an increase in transcript abundance. These data demonstrate that even with the potential for polyploidy, relative quantification is applicable for measuring gene expression changes in extremophiles. Representative figure of data used to obtain transcript-to-gene ratios for gyrB, polA2, and rpoC. Samples were collected at multiple time points from anaerobic cultures supplemented with glucose as the carbon source and nitrate as the terminal electron acceptor. Transcript and gene copies per mL of sample were calculated via absolute quantification using external calibration curves. Error bars represent the standard deviation derived from triplicate qPCRs from three biological replicates. doi:10.1371/journal.pone.0131015.g002

Absolute Quantification to Identify Specific Fluctuations
The data obtained via absolute quantification was then used to identify the specific conditions and times under which gene expression fluctuated. Gene copy number was similar among the three genes under each condition (representative image of transcript and gene copy numbers per mL for each of the three genes is shown in Fig 2, transcript and gene copy numbers per mL for all three genes under all conditions shown in S5-S13 Figs), corroborating genome annotation analysis that showed each gene existing as a single copy within the genome. PCR efficiencies (PCR E) were similar among the three assays (Table 2) and so should result in similar gene copy numbers under a specific treatment barring unusual changes in cell physiology. Therefore, ratio fluctuations could be attributed to changes in transcript copy numbers.
There was no significant difference among gyrB ratios at either 26 or 45 h of growth under aerobic or anaerobic cultures, nor was there any significant difference in ratios across time under any of these three conditions (Fig 3A-3C). These results corroborated with those obtained with BestKeeper based solely on C T values. There was also no significant difference between gyrB ratios when cells were grown anaerobically with iron or nitrate as the terminal electron acceptor (Fig 3D). In comparing the ratios of gyrB for cells grown in an oven compared to a microwave, there was no significant difference at 8 h, but there was at 14 h (t-test, p = 0.00) (Fig 3E and 3F); additionally, there was no significant difference in gyrB ratios over time during microwave incubation (one-way ANOVA, p = 0.06), but there was during oven incubation (one-way ANOVA, p = 0.00).
Ratios for rpoC were also found to differ significantly at multiple time points when comparing aerobic (TYG) and anaerobic (with glucose or lactate as carbon source in a minimal medium) growth (26 h, one-way ANOVA, p = 0.03, 45 h one-way ANOVA, p < 0.00) ( Fig  4A-4C). Additionally, there was not a significant difference in rpoC ratios between cells grown with different terminal electron acceptors (Fig 4D). Ratios for rpoC followed a similar trend to that of gyrB for oven and microwave growth, with a significant difference at 14 h (t-test, p = 0.00) but not 8 h (Fig 4E and 4F). There was also a significant difference in rpoC ratios across time when cells were grown in an oven (one-way ANOVA, p = 0.02) or with microwave heating (Kruskall-Wallace, p = 0.04).
PolA2 differed significantly at both 26 h and 45 h when comparing aerobic and anaerobic growth (one-way ANOVA, p = 0.01, p = 0.00, respectively), and also across time when grown anaerobically in minimal medium with glucose as the carbon source (one-way ANOVA, p = 0.01) (Fig 5A-5C). However, there was no significant difference across time when cells were grown anaerobically with lactate as the sole carbon source, or aerobically in TYG (Fig 5B and  5C). There was also no significant difference in the ratios when cells were grown anaerobically with different electron acceptors ( Fig 5D). As with gyrB and rpoC, polA2 ratios were not significantly different after 8 h growth in oven or MW, yet were after 14 h (t-test, p = 0.00). PolA2 ratios were also found to differ significantly across time under both oven and microwave incubation (one-way ANOVA, p = 0.01; Kruskall-Wallace, p = 0.04, respectively) (Fig 5E and 5F).

Conclusions
In this study, we designed and validated assays for a suite of reference genes for use in T. scotoductus expression studies under multiple culture conditions. Based on their unique biomolecular properties that have allowed them to adapt to extreme environments, the question arises as to whether mesophilic reference genes are suitable for use with extremophiles. Under aerobic growth and anaerobic growth with different electron donors and nitrate as the terminal electron acceptor, 19900, polA2, gyrA, and gyrB ranked as the most stable genes. Under conditions of anaerobic growth with varied electron acceptors, 19900, infA, pfk, gyrA and gyrB were the most stable. Three genes displayed a SD < 1.00 under aerobic growth using different heating methods: gyrA, gap, and gyrB. Overall, the results demonstrate that when assimilating all culture conditions, gap, gyrA, and gyrB are the most stable (SD < 1.00) and thus the most appropriate if using a combination of conditions. The consistent stability of gyrA is similar to the results obtained in reference gene studies with mesophiles. We also demonstrated that rpoC does not serve as a reliable reference gene in thermophiles, due to its expression instability across all culture conditions. Absolute quantification with three of the genes demonstrated that even with the potential for polyploidy, relative quantification is applicable for use in RT-qPCR expression studies with extremophiles.
Supporting Information S1 Data. Raw C T values and intra-assay variations for 11 reference genes assessed in this study. For genes (gyrB, polA2, rpoC) assessed via absolute and relative quantification, "DNA" refers to qPCRs in which DNA served as template and was used to calculate gene copies per mL; "RNA" refers to qPCRs in which cDNA served as template and was used to calculate transcript copies per mL. (XLSX) Fig 5. Correlation between polA2 transcript C T values and gene expression as derived from the transcript-to-gene copy number ratios. A standard curve was constructed and used to calculate the transcript copies per mL of sample and gene copies per mL of sample ("transcript-to-gene ratio") under the different culture conditions described in the text. Ratio and C T values obtained for polA2 under: anaerobic growth with (A) glucose or (B) lactate as carbon source compared to (C) small-volume aerobic growth in TYG; (D) anaerobic growth with glucose as carbon source and nitrate or iron as terminal electron acceptor; and aerobic growth in TYG using (E) oven or (F) microwave heating. Glu = glucose, Lac = lactate, MW = microwave heating. a = indicates statistically significant difference (p < 0.05) among different culture conditions at specific time points; b = indicates statistically significant difference (p < 0.05) across time under one culture condition. A 700 base pair region of the gyrase subunit B gene from the Thermus scotoductus genome was amplified using gene-specific primers. The resulting PCR product was purified and used in the construction of an external calibration curve. The number of copies per μL was calculated using the equation: X g μL -1 DNA/(PCR amplicon length x 660) x 6.022 x 10 23 . Serial 10-fold dilutions were prepared that spanned eight orders of magnitude ranging from 1 x 10 2 to 1 x 10 9 copies. The standard curve was determined by plotting the log of the calculated copy number against the cycle at which fluorescence for that sample crossed the threshold cycle. The slope, y-intercept, r 2 value, and PCR efficiency of the gyrB assay is provided in S1 Table. (TIF) T. scotoductus was grown anaerobically in basal medium with glucose as the carbon source and iron (Fe) or nitrate (NO 3 ) as the terminal electron acceptor. Transcript and gene copies per mL of sample were calculated via absolute quantification using external calibration curves as described in the text. Error bars represent the standard deviation of triplicate qPCRs from three biological replicates.
(TIF) S11 Fig. polA2 transcript and gene copies over time comparing anaerobic versus aerobic culture conditions. T. scotoductus cultures were grown anaerobically in a basal medium with lactate as the carbon source and nitrate as the terminal electron acceptor (top panel) or aerobically in TYG medium (bottom panel). Transcript and gene copies per mL of sample were calculated via absolute quantification using external calibration curves as described in the text. Error bars represent the standard deviation of triplicate qPCRs from three biological replicates. (TIF) S12 Fig. polA2 transcript and gene copies over time comparing oven versus microwave heating. T. scotoductus cultures were grown aerobically in TYG using oven (top panel) or microwave (bottom panel) heating. Transcript and gene copies per mL of sample were calculated via absolute quantification using external calibration curves as described in the text. Error bars represent the standard deviation of triplicate qPCRs from three biological replicates. (TIF) S13 Fig. polA2 transcript and gene copies after 24 h comparing terminal electron acceptors. T. scotoductus was grown anaerobically in basal medium with glucose as the carbon source and iron (Fe) or nitrate (NO 3 ) as the terminal electron acceptor. Transcript and gene copies per mL of sample were calculated via absolute quantification using external calibration curves as described in the text. Error bars represent the standard deviation of triplicate qPCRs from three biological replicates. (TIF) S1 Table. Standard curve details for gyrB, polA2, and rpoC assays, including primers used in standard construction, slope, y-intercept, r 2 , and PCR efficiency.