Selection and Evaluation of Reference Genes for Expression Studies with Quantitative PCR in the Model Fungus Neurospora crassa under Different Environmental Conditions in Continuous Culture

Neurospora crassa has served as a model organism for studying circadian pathways and more recently has gained attention in the biofuel industry due to its enhanced capacity for cellulase production. However, in order to optimize N. crassa for biotechnological applications, metabolic pathways during growth under different environmental conditions must be addressed. Reverse-transcription quantitative PCR (RT-qPCR) is a technique that provides a high-throughput platform from which to measure the expression of a large set of genes over time. The selection of a suitable reference gene is critical for gene expression studies using relative quantification, as this strategy is based on normalization of target gene expression to a reference gene whose expression is stable under the experimental conditions. This study evaluated twelve candidate reference genes for use with N. crassa when grown in continuous culture bioreactors under different light and temperature conditions. Based on combined stability values from NormFinder and Best Keeper software packages, the following are the most appropriate reference genes under conditions of: (1) light/dark cycling: btl, asl, and vma1; (2) all-dark growth: btl, tbp, vma1, and vma2; (3) temperature flux: btl, vma1, act, and asl; (4) all conditions combined: vma1, vma2, tbp, and btl. Since N. crassa exists as different cell types (uni- or multi-nucleated), expression changes in a subset of the candidate genes was further assessed using absolute quantification. A strong negative correlation was found to exist between ratio and threshold cycle (CT) values, demonstrating that CT changes serve as a reliable reflection of transcript, and not gene copy number, fluctuations. The results of this study identified genes that are appropriate for use as reference genes in RT-qPCR studies with N. crassa and demonstrated that even with the presence of different cell types, relative quantification is an acceptable method for measuring gene expression changes during growth in bioreactors.


Introduction
Neurospora crassa is a filamentous fungus that has served for nearly 50 years as the model organism for studying circadian pathways and the molecular clock in eukaryotes [1].It is classified under the phylum Ascomycota and is widely distributed in nature [1], typically found growing on vegetation killed by fire [2].As such, it is able to synthesize and secrete high levels of all main enzyme types involved in lignocellulosic biomass degradation and is also able to convert pentose sugars, cellulose polymers, and agro-industrial residues to ethanol [3].Its genome is predicted to contain twice as many cellulase genes in comparison the current state-of-the-art industrial strain used in biofuel production, Trichoderma reesei [3].Additionally, recent work has shown that the variety of molecular, genetic, and biochemical techniques available for N. crassa can advance the analyses of fungal deconstruction of plant biomass [4].Therefore, great interest in the use of N. crassa for biotechnological applications has sparked in recent years.
Gene expression profiling is an effective means of studying an organism's response to its environment, a necessary component in identifying the genes and regulatory mechanisms associated with biological processes in N. crassa such as circadian pathways and important biotechnological applications such as cellulase degradation.Recent work in our lab has utilized N. crassa maintained in continuous, long-term bioreactor cultures in an effort to understand the triggers of cellulase production.Most N. crassa studies are based on either racetube or batch culture conditions.Understanding the complexity of N. crassa growth in bioreactors is lacking and more importantly, the identification of optimal reference genes for molecular analysis needs to be addressed.
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.With this technique, two quantification strategies are available to measure the level of expressed genes: absolute and relative quantification.Absolute quantification relates the fluorescent PCR signal of the sample to input copy number using an external calibration curve, while relative quantification measures the relative change in mRNA expression levels [5].Relative expression is based on normalization of target gene expression to an internal control gene, often referred to as a ''reference'' or ''housekeeping'' gene, that should not vary in the cells under investigation or in response to environmental or experimental conditions [6].The selection of a suitable reference gene is critical for gene expression studies using relative quantification [5].However, these reference genes, thought to be constitutively expressed and minimally regulated, are commonly used for expression profiling in RT-qPCR assays without a priori conformation that mRNA levels do not change under the experimental conditions being investigated.The use of reference genes necessitates that they be validated for the specific experimental set-up and it is typically necessary to choose more than one [7].Work by Vandesompele et al [6] demonstrated that the typical procedure in which a single reference gene is used for normalization can lead to relatively large errors [6], and it is now common to employ a combination of genes.Thus, in order to understand the metabolic processes when N. crassa is grown under different environmental conditions, the selection and evaluation of reference genes must be addressed.
N. crassa is a morphologically complex multicellular organism, possessing at least 28 distinct cell types [8] between the vegetative (asexual) and sexual phases.In bioreactors, it is maintained in the vegetative phase.The primary cells occurring in the vegetative phase include hyphae, macro-and microconidia, all of which are comprised of several different types [8].One of the biological differences among these cell types is the number of nuclei per cell: the macroconidia are multinucleated, possessing on average between 3-6 nuclei per cell, while microconidia possess only a single nuclei (Figure 1).This life cycle trait must be taken into account when performing gene expression studies, particularly under continuous, long-term culture conditions such as those conducted in bioreactors.In general, our bioreactors are dominated by clusters of aggregated colonies dispersed throughout the vessel, similar in cell type composition to the mats observed in batch cultures yet comprised of more discrete units.Therefore, it is necessary to determine whether changes in gene expression are the result of true fluctuations in transcripts, or whether an increase or decrease in transcripts (i.e. gene expression) is actually a reflection in gene copy number fluctuation, which would occur when cells are a mix of macro-and micro-conidia.One means of assessing gene expression changes is to determine the number of transcripts normalized to a certain value such as cell number, volume of sample, or mass.This number can then be related to the gene copy number normalized to the same units to provide a ratio of the transcript to gene ratio.These data standardize the transcripts to the gene copy number and so determine whether the gene is changing expression, or whether transcript variation is a result of gene variation.This type of situation applies to cultures that occur as multiple cell types during which some morphologies are multinucleated.
While numerous gene expression studies have been performed with N. crassa, little to no information is provided as to reference gene validation, and the extent of gene expression fluctuation as a result of changes in gene copy number (i.e. as cells progress from hyphae to macroconidia) has not been addressed.Additionally, genes that are commonly used as reference genes in other organisms have been shown to be under clock control in N. crassa.In general, the actin gene is frequently used as a normalizer in gene expression studies [9][10][11], though in most cases there is no mention of stability verification.Therefore, the overall objective of this study was to examine a suite of genes for the potential to serve as appropriate reference genes under a variety of environmental conditions during growth in N. crassa.We evaluated a set of 12 genes under different environmental conditions using a combination of both relative and absolute quantification.Using this combination allowed for a more comprehensive assessment of fluctuations in expression, as the transcript copy number per mL was calculated and its ratio to gene copy number determined.Additionally, the well-established software programs BestKeeper [12] and NormFinder [13] were used to determine expression stability and rank of all genes to assess their suitability as reference genes.

Organism Maintenance
N. crassa cultures were maintained on agar slants (1.5% w/v sucrose with Vogel's medium [14]) at 220 ˚C for up to 1 month.Slants were prepared following standard protocols.Briefly, the medium comprised of Vogel's salts, 1.5% sucrose, and 1.5% agar was autoclaved and cooled to #50 ˚C whereupon the trace minerals and biotin solutions (1.0610 25 M) were added.Agar was distributed among test tubes and cooled at an angle, resulting in a ''slanted'' agar medium.Strain X-661 (ras-1bd mutant with the frq:luc transcriptional reporter) was used for these studies.

Bioreactor culture conditions and cell harvesting
New slants were inoculated from the frozen stock culture at the start of each experiment.N. crassa X-661 grew on slants for 2-3 days at 30 ˚C in constant darkness.Conidia were isolated from slants by washing the slant culture with ,1 mL of Vogel's medium [14] lacking a carbon source.The collected material was centrifuged (10,000 rpm, 1 min) and the resulting conidial pellet was washed 3 times with Vogel's medium (no C-source).The final pellet was re-suspended and inoculated into 60 mL of Vogel's medium with 1% glucose, 0.2% Junlon, (100,000 MW, Sigma Aldrich) and 0.008% Antifoam 204 (Sigma Aldrich).The culture was grown at 30 ˚C with shaking for 24 hrs with light.This culture was inoculated into a BioFlo/CelliGen 115, 1.3 L bioreactor (Figure 2) containing the same medium composition and grown in batch for 24 hrs with light.Continuous flow into and export from the culture was started as the cells reached the exponential growth phase.When the cells reached steady state (,11 hrs later) the light remained on for 11 hrs and then was shut off.At this point, the cells were grown in one of three conditions: continuous all-dark conditions, under light/ dark cycling, or under temperature flux (25 ˚C/37 ˚C) in constant darkness.The light/dark and temperature flux oscillated with 22-hr periods (11 hrs light, 11 hrs dark) in keeping with the period length of N. crassa circadian oscillations.During the experiments, N. crassa cultures were maintained in steady-state conditions, with pH 5.5 and 10% dissolved O 2 .The O 2 was bubbled through the medium at 1.0 L/min.A dynamic agitation program was used to keep the cells suspended in bioreactor cultures.The program controlled the rotation rate of the propeller in the reactor and was set up as a 10-min cycle using the BioCommand Batch Control Plus BioProcessing software (New Brunswick Scientific, #M1326-0051, revision B).The propeller rotated at 1,000 rpm for 1 min, 100 rpm for 1 min, and then 400 rpm for 8 min.This program was started ,1-2 hrs following bioreactor inoculation.The temperature was either maintained at 25 ˚C or cycled (with 11-hr periods) between 25 ˚C and 37 ˚C.The bioreactors were equipped with an external heating jacket and an internal cooling finger and temperature probe to accurately maintain the temperature.A temperature cycling program (BioCommand software) was used for varied temperature experiments.The growth state of the cells (i.e., log phase, steady-state) was monitored using the CO 2 output.A CO 2 gas sensor (Vernier) was used to detect CO 2 produced by the culture.The CO 2 data was recorded using the LoggerLite software (version 1.6.1).Samples were removed under darkroom light conditions (safelight) and the 1.5 mL aliquots were flash frozen in centrifuge tubes and stored at 280 ˚C.Biological replicates, in the form of individual bioreactors, were performed for each type of environmental condition, with a total of 8-10 time points collected per condition.

DNA extraction
Genomic DNA was extracted from biological replicates using the DNeasy Plant Mini kit (Qiagen, Valencia, CA) following the manufacturer's instructions, with slight modifications.Briefly, cells from 1 mL of a bioreactor sample were disrupted using a mortar and pestle under liquid nitrogen.Following addition of 400 mL Buffer AP1 and 4 mL RNase A, the slurry was ground an additional 1-2 min.The entire volume was then transferred to a 1.5 mL microfuge tube, vortexed vigorously, and 400 mL of this slurry removed for DNA extraction.The remaining steps followed the manufacturer's protocol.DNA was eluted in nuclease-free water rather than Buffer AE.DNA concentration and purity were assessed using the Nanodrop-2000 (Thermo Fisher Scientific, Wilmington, DE).DNA used in quantitative PCR (qPCR) assays was diluted 1:5 in nuclease-free water prior to use.

RNA isolation and cDNA synthesis
Total RNA was extracted frozen from a 1 mL bioreactor sample using the RNeasy Plant Mini kit (Qiagen, Valencia, CA) following the manufacturer's protocol, with slight modifications.Briefly, cells were disrupted using a mortar and pestle under liquid nitrogen.Four hundred fifty microliters of buffer RLC was added and the slurry mixed until pipetteable.The sample was vortexed vigorously and 450 mL added to the QIAshredder spin column.Subsequent steps followed the manufacturer's protocol.The DNase digestion was extended to 1 hr.An additional spin at top speed (16,000 rpm) for 1.5 min was included immediately prior to elution in 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).The protocol was optimized to include less template RNA and a larger overall reaction volume, as preliminary screening with several of the potential reference genes revealed greater fluctuations in C T s when using the recommended 2 mg in a 20 mL volume (data not shown).The following protocol was used for all reverse transcriptase reactions: total RNA was first diluted 1:5 in nuclease-free water, and 500 ng added to the reaction.Each reaction contained: 2 mL 20x enzyme mix, 2x RT buffer, a volume equivalent to 500 ng total RNA, and brought to a final volume of 40 mL with nuclease-free water.Controls consisted of all components except the enzyme, in which case an equal volume of nuclease-free water was added.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.Prior to use in qPCR, samples were again diluted 1:5 in nuclease-free water.

Candidate reference gene primer design and validation
Quantitative PCR assays were developed for a suite of candidate reference genes (Table 1) based on the SYBR Green chemistry [15].For each assay, primers were designed that targeted an ca.100-130 bp region of the coding sequence of each gene (Table 1).When possible, primers were designed to span exon-exon boundaries.Both coding and genomic sequences of each gene were obtained from the N. crassa database located on the Broad Institute website (http://www.broadinstitute.org/annotation/genome/neurospora/GeneIndex.html).Primers were designed using the Primer3 software (http://bioinfo.ut.ee/primer3/), with an annealing temperature of 60 ˚C.Primers were examined for potential secondary structures using Mfold [16].The potential for self and cross dimer formation was examined for each of the primer sets using an online oligonucleotide analysis tool (https://www.operon.com/oligos/toolkit.php).Self-complementarity was determined using the 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 was determined using LinRegPCR [17].

Construction of external standards for absolute quantification
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 though the creation of external calibration curves.Plasmidbased standards were constructed for the actin (act), ATP synthase lyase c subumit (asl), adenosine kinase (adk), and vacuolar ATPase subunit 3 (vma3) genes and used for absolute quantification of both transcript and gene copy number.Standards were created via amplification from genomic DNA using genespecific primers (Table S1).DNA was used rather than RNA as it has been shown that external calibration curves constructed with DNA rather than RNA serve as Sequence identity was confirmed using the standard nucleotide BLAST function on the NCBI website, with organism specificity set to ''Neurospora.''The number of copies per mL was calculated using the equation: X g mL 21 DNA/[PCR amplicon + plasmid length6660])66.022610 23.The plasmid stock solution was diluted in nuclease-free water to a working concentration of 5610 8 copies per mL.Serial 10-fold dilutions were used to create external calibration curves that spanned eight orders of magnitude ranging from 5610 7 to 5610 21 copies per mL.Two mL of each dilution were used per reaction, yielding a standard curve that ranged from 1610 1 to 1610 8 copies.The standard curve for each assay was obtained by plotting the log of the calculated copy number against the cycle at which fluorescence for that sample crossed the threshold value (cycle threshold, C T ) (Figures S1, S2, S3, and S4).The slope, y-intercept, r 2 value, and PCR efficiency of each assay can be found in the Supporting Information (Figures S1, S2, S3, and S4).

Quantitative PCR
The Fast SYBR Green Master Mix (Applied Biosystems, Foster City, CA) was used for all qPCR assays.Each reaction contained: 10 mL 2x Fast SYBR Green master mix, 150-300 nM each forward and reverse primer, 5 ng (2 mL) template cDNA, and brought to a final volume of 20 mL with nuclease-free water.Reactions were conducted in 96-well plates on a ViiA 7 with the Fast 96-well block system (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 96-well reaction plates (0.1 mL with barcode) sealed with MicroAmp optical adhesive film.Notemplate 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 C T values and intra-assay variations for both relative and absolute quantification assays are provided in Data S1 and Data S2.

Data analysis
Three types of environmental conditions were used to determine the appropriateness of a suite of genes for use as reference genes: total darkness, in which samples were collected over a time span of 70 hrs; light/dark cycling, in which samples were collected from alternating periods of light and dark over the span of ca.70 hrs; and temperature fluctuations in darkness, in which samples were collected from alternating cycles of 25 ˚C and 37 ˚C over the span of ca.70 hrs.Expression levels of all twelve candidate reference genes were initially determined from their C T values.The expression stability of all candidate reference genes under each of the different conditions and among all conditions combined was assessed using the BestKeeper v. 1 and Normfinder software programs.For data input into NormFinder, C T s of each gene at each time point were converted to relative quantities using the equation Q5E DC T , where DC T 5C T highest abundant sample -C T sample, E5PCR efficiency, and Q5relative quantity.
The number of transcripts and gene copies per mL were calculated via external calibration curves for act, adk, asl, and vma3.Transcript or gene copy numbers from total RNA or DNA extracted from 1 mL bioreactor samples were calculated automatically from the regression lines by the ViiA 7 System software.Gene copies were calculated using an equation described previously [19], in which gene copies per mL sample5(gene copies per reaction mix x volume of DNA [mL])/(DNA per reaction mix [mL]6mL sample used) and adjusted to account for the 1:5 initial dilution factor.The following equation was used to calculate transcript copies per mL: [(copies per mL per reaction mix6RNA dilution factor)/(volume mL total RNA added to RT reaction/RT reaction volume mL6cDNA dilution factor]6vol total RNA eluted (mL) (Data S2 and Data S3).The transcript to gene ratio was determined by dividing transcript copies per mL by gene copies per mL.The Shapiro-Wilkes test was used to test for normality of both the ratio and C T values (Table S2).If both data sets passed the test for normality, the Pearson product moment correlation test was applied to examine the strength of the correlation between the two.If one of the data sets failed the test for normality, the Spearman rank order correlation test was applied.

Results and Discussion
Quantitative PCR assay design, optimization, and general characteristics Typical reference genes include GAPDH, albumin, actins, tubulins, and small and large subunit rRNA [5], as they are present in all nucleated cell types and are necessary for basic cell survival.However, numerous studies, across different organisms, have shown that the typical reference genes are regulated and vary under different conditions [5].Based on the experimental evidence thus far, it appears that no reference gene retains sufficient overall expression stability suitable for all contexts and that the best candidates differ between tissue type [20] and environmental condition.Therefore, it is recommended that the individual investigator choose a reference gene that is most stably expressed under their specific conditions, such as tissue type or stage of development, and should be impervious to perturbations such as external stimuli, disease, or genetic modification [5,[20][21][22].
Several approaches exist for qPCR data normalization.One approach recently implemented has been to select genes from a genome-wide background derived from large sets of microarray data.This approach has been used for Arabadopsis, multiple human cell lines and cancers, and mouse [20].However, while N. crassa serves as a model organism for studying eukaryotic circadian pathways, it is not included in this database.Additionally, studies utilizing whole-transcriptome profiling such as microarrays and RNA sequencing are able to select appropriate reference genes based on the transcriptomic profiling.For example, genes encoding an adenosine kinase and the ATP synthase gamma chain have been used as reference genes in the filamentous fungus Penicillium glabrum under conditions of heat stress, as their expression was shown via initial microarray analysis to be unaffected by this condition [23].Additionally, a recent study evaluated 12 potential reference genes in the yeast S. cerevisiae based on pooled data from publicly-available transcriptome studies.The majority of selected genes were not previously reported as reference genes and the results revealed that these newlyidentified genes outperformed the commonly-used reference genes [24].There is extremely limited gene expression data with N. crassa growth in bioreactors, such that large data sets from which to select potential housekeeping genes do not exist.Also lacking are publicly-available whole-transcriptome data sets such as has been used for other organisms when selecting reference genes.Therefore, potential reference genes were selected based on the following features: previous transcriptomic-based studies from N. crassa (acl) [25] or other fungi (adk, asl) [23] in which the expression of these genes remained unchanged; via selection from the Broad Institute N. crassa database (http://www.broadinstitute.org/annotation/genome/neurospora/MultiHome.html) with genes for which no expression data existed (tbp), or from previous studies of genes whose results pertaining to their function and regulation indicated the potential for use as a reference gene (vma1, vma2, vma3, trk) [26,27] Additionally, N. crassa genes employed as reference genes in previous studies were also evaluated here (act, btl, l6, frh) [28], albeit with newly-designed primers for use in the btl and frh assays.
The compliance of the RT-qPCR assays with the MIQE (Minimum Information for Publication of Quantitative Real-Time PCR Experiments [29]) guidelines is shown in the MIQE checklist (Table S3).All assays developed in this study were specific for their intended gene target, as indicated by a single peak in the melt curve analysis (Figure S5).Additionally, all possessed similar PCR efficiencies (Table 1), so that if all had displayed equal stability, they could have been used interchangeably in relative quantification methods.A preliminary assessment of the 12 reference genes showed that average C T values, when assimilating the data from all genes, ranged from 21.47-30.34under light/dark fluctuations (Figure 3A), 22.53-30.72under all dark conditions (Figure 3B), and 19.93-28.73under temperature flux (Figure 3C).Based on this initial analysis, the two genes displaying the highest average C T values, frh and trk, were discarded from further analysis.This was due to the recommendation of reference genes being expressed at comparable levels to both other reference genes and target genes [30].Additionally, frh has been shown to be under clock control [31].

Reference gene stability under light/dark conditions
The expression stability of candidate reference genes was initially estimated 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.With this program, reference genes can be ordered from the most stably expressed, exhibiting the lowest variation, to the least stable.In addition to stability (SD and CV) values, additional key parameters include the Pearson correlation coefficient (r), in which pair-wise correlation analyses are performed to estimate inter-gene relationships among all reference genes, and their associated probability (p) values.Thus, r and p describe the relation between the BestKeeper index and each gene.The candidate reference genes were then ranked as to their suitability under one environmental condition and across all three using the program NormFinder.This program utilizes a model-based approach that, in addition to assessing overall expression variation, allows for the samples to be assigned to different groups, allowing for the intra-and intergroup expression variation of the candidate genes to be assessed; this grouping feature was used when assessing reference gene stability under light/dark and temperature cycling.It was not employed for all-dark growth at constant temperature.The comprehensive results of the BestKeeper and NormFinder programs are provided in Data S4 and Data S5, respectively.
The greatest variability among the three environmental conditions was observed under light/dark fluctuation, as evidenced by both the range in C T s and the SD values obtained with BestKeeper (Table 2).In general, most genes displayed overall stability with respect to light/dark cycling.Three genes (btl, asl, and vma1) displayed the greatest expression stability (SD51.10,1.13, and 1.21, respectively), while l6 and acl were the least stable (SD51.41).However, based on the Pearson correlation coefficients (r), they correlated well with the other genes and thus were retained for further analysis per previous methods [12].Overall, genes displayed similar Pearson correlation coefficients, ranging between 0.936-0.997.NormFinder identified vma1 as the most stable gene (Table 2), with the best combination of two genes being vma1 and asl.The similarity in stability values between vma1 and asl was also reflected in the intragroup variation, which was 0.001 in the light for both.However, vma1 expression was less influenced under dark conditions than asl, with an intragroup variation of 0.002 versus the 0.016 calculated for asl (Table 2).

Reference gene stability under all-dark conditions
Based on the SD as calculated by BestKeeper, all genes were stable during complete darkness, as the SD for all genes was ,1 (Table 3).The only exception to this was vma3, in which the SD51.07;however, the Pearson correlation coefficient was acceptable, and in fact higher than some of the other genes that displayed a lower SD than vma3.Genes displaying the lowest SD, and thus the greatest stability, were vma1, btl, and tbp.All genes were significantly correlated to the BestKeeper index (p,0.001),and the Pearson correlation coefficient was very similar among genes.In general, the analysis showed a strong correlation for all 10 genes.NormFinder identified tbp, vma2, and vma1 as being the most stable of the genes under conditions of complete darkness, with tbp and vma2 recording almost identical stability values (0.141 and 0.143, respectively) (Table 3).

Reference gene stability during temperature fluctuation
All genes were stably expressed under conditions of temperature fluctuations, as reflected by the SD ranging from 0.501-0.865.Overall, SDs ranged between 0.501-0.741,with the exceptions of l6 (0.796) and acl (0.865).Based on the SD, asl (0.501), btl (0.616), and vma1 (0.622) exhibited the greatest stability (Table 4).As had been observed for all-dark conditions, all genes correlated very strongly with the index, with most displaying r.0.9.Under conditions of fluctuating temperature, Normfinder identified vma1 as the most stable, with the best combination being vma1 and act (Table 4).vma1 displayed extremely low variation under 25 ˚C (0.0001) with an increase in variance observed at 37 ˚C for both vma1 and act (0.026 and 0.018, respectively).

Reference gene stability across all conditions
Assimilating the results of both BestKeeper and NormFinder under each individual condition, vma1 always ranked among the top three genes in terms of stability; asl and btl also consistently ranked in the top three in at least one of the analyses.Also, in keeping with their biological role as subunits of the same enzyme, vma1, vma2, and, under several conditions, vma3 ranked closely together in terms of stability, suggesting they can be used interchangeably.However, in order to accurately identify the most stable genes under all three conditions, all data were combined and analyzed with both BestKeeper and Normfinder.The results of BestKeeper revealed that 4 of the 10 genes displayed acceptable stability across all conditions (reflected by a SD,1): btl, tbp, vma2, and vma1 (Table 5).
Conversely, the two genes frequently used for normalization, l6 and act, displayed the greatest variance (Table 5).The most appropriate gene for use across all conditions as identified via NormFinder was vma1, with the best combination consisting of vma1 and asl, with a combined stability value of 0.089.tbp and vma2 were the next most-stable genes among the 10, in agreement with the results obtained with BestKeeper (Table 5).

Absolute quantification of transcript and gene copy numbers
In bioreactors, N. crassa exists as different vegetative cell types, some of which are multi-nucleated.Based on the parameters [6,12] defined for appropriate reference gene selection, all 10 candidate reference genes displayed fluctuations in expression.This prompted the question: Were fluctuating C T s due to changes in gene expression, or were gene (DNA) copy numbers changing, so that if transcript levels remained relatively stable, it was manifested as expression fluctuations?Therefore, expression changes in a subset of the candidate genes was assessed using absolute quantification.The transcript (RNA) and gene (DNA) copy numbers of act, adk, asl, and vma3, were calculated and used as an additional means of examining gene expression.These four were selected based on combination of factors since they are involved in different biological processes within the cell and ranged in their stability values under the conditions examined.For example, asl frequently ranked as one of the most stable, while adk was one of the least stable.External 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 over time (representative image in Figure 4, plots of transcript and gene copy numbers for all genes under all conditions Figures S6, S7, S8, and S9).This ratio indicated whether C T fluctuations were due to changes in transcript abundance, indicative of changes in gene expression, or whether it was a result of changes in gene copy number, such as could occur in multinucleated cells.
Plotting both the RNA: DNA ratio and C T values over time indicate that C T changes are due to fluctuations in transcript copy number (Figure 5).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 (i.e. as the ratio increases, the C T decreases).With all four genes under all three conditions, a strong negative correlation was found to exist between ratio and C T values, demonstrating that C T changes serve as a reliable reflection of transcript, and not gene copy number fluctuations (Table 6).For example, under conditions of light/dark cycling, the act C T value decreased at the second time point (Figure 5A), indicative of an increase in transcripts.If this were due to changes in gene (DNA) copy number, the RNA: DNA ratio would not increase.However, an increase in ratio indicates an increase in transcript abundance.As statistically significant negative correlations were found to exist between the transcript: gene ratio and C T values among the four genes under all conditions, these data demonstrate that even with the presence of different cell types, relative quantification is an acceptable method for measuring gene expression changes in N. crassa during growth in bioreactors.

Conclusions
The work presented here mirrors the results obtained in recent studies with multiple organisms in which the most appropriate reference genes for use in RT-qPCR were not those commonly employed when using the relative quantification method, underscoring the importance of proper reference gene evaluation and selection.The increased use of RT-qPCR coupled with the realization of the importance of selecting the appropriate reference gene has led to dedicated studies on reference gene selection in both model and non-model organisms, including rats [32,33], zooplankton [34], fruit crops [35] (Litchi chinensis), and insects [36][37][38][39].In some cases, no single gene was found to be stable under all conditions, with stability affected by both endogenous and exogenous factors such as cell type, developmental stage, diet, and environmental stressors [38,39].Additionally, reference gene studies have been undertaken in recent years with multiple genera of filamentous fungi.A study evaluating potential reference genes in Trichoderma reesei, similar in industrial importance to N. crassa due to its high secretory capacity for hydrolases, identified sar1, encoding a small GTPase, as most stable, whereas the gene coding for actin (act), did not rank among the best validated ones [40].In the filamentous fungus Fusarium graminearum, an economically-  important plant pathogen, evaluation of 15 genes previously identified as housekeeping genes or those selected from the whole transcriptome sequencing data under different culture conditions found the most appropriate reference genes to vary based on growth stage and/or toxin production [41] and were comprised of genes coding for a cell surface flocculin, ubiquitin C-terminal hydrolase, the eukaryotic translation initiation factor 1 alpha, and the mitochondrial ribosomal protein S16.
In this study, we identified and validated new reference genes for use in N. crassa expression studies under multiple environmental conditions, including long-term growth in constant darkness, light/dark cycling, and temperature flux in constant darkness.In addition to different environmental conditions, we assessed the time-course variation in reference gene expression by sampling at frequent intervals over a period of .60 hours.These diverse conditions and continuous culture conditions span primary developmental (mycelia, macro-and microconidia) and/or physiological stages, thus increasing the reliability of the reference genes selected in this study.Overall, the results demonstrate that, under the culture conditions described here, vma1, vma2, tbp, and btl serve as the most appropriate reference genes for use in qPCR.We also demonstrated that l6, and to a lesser extent act, often used in the relative quantification method, are less suitable as reference genes in comparison to multiple other genes.fold dilutions were used to create external calibration curves that spanned eight orders of magnitude ranging from 1610 1 to 1610 8 copies.The standard curve for each assay was obtained by plotting the log of the calculated copy number against the cycle at which fluorescence for that sample crossed the threshold value.The slope, yintercept, r 2 value, and PCR efficiency of the actin assay is provided in the figure.

Figure 1 .
Figure 1.Generalized lifestyle of N. crassa during which different cell types occur.Center image illustrates bioreactor set-up.Within the bioreactors utilized in this study, N. crassa was present as both multiand uninucleated cells.doi:10.1371/journal.pone.0112706.g001

Figure 3 .
Figure 3. Boxplot of C T values of candidate reference genes during light/dark cycling, temperature fluctuation, and continuous darkness.A line across the box depicts the median.The box indicates the 25th and 75th percentiles.Whiskers represent the 90th and 10th percentiles.doi:10.1371/journal.pone.0112706.g003

Figure 4 .
Figure 4. Representative figure of data used to obtain transcript-to-gene ratios.act all-dark 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 biological replicates.doi:10.1371/journal.pone.0112706.g004

Figure 5 .
Figure 5. Correlation between gene expression and C T values as derived from absolute quantification of transcript and gene copy numbers.Standard curves were created for the act, adk, asl, and vma3 genes and used to calculate the transcript-to-gene ratios under conditions of light/dark cycling, temperature fluctuation, and continuous darkness.Top panel depicts light/dark cycling, with yellow highlight indicating light on.Bottom panel displays conditions of temperature flux between 25 ˚C and 37 ˚C, with 37 ˚C illustrated by orange highlight.doi:10.1371/journal.pone.0112706.g005 Figure S2.Standard curve of adk assay, constructed as described for the act assay.doi:10.1371/journal.pone.0112706.s002(TIF) Figure S3.Standard curve of asl assay, constructed as described for the act assay.doi:10.1371/journal.pone.0112706.s003(TIF) Figure S4.Standard curve of vma3 assay, constructed as described for the act assay.doi:10.1371/journal.pone.0112706.s004(TIF) Figure S5.Melt curve profiles of the 12 reference gene assays evaluated in this study.doi:10.1371/journal.pone.0112706.s005(TIF) Figure S6.act transcript and gene copies over time under conditions of light/ dark cycling, temperature flux, and continuous darkness.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 derived from triplicate qPCRs from biological replicates.doi:10.1371/journal.pone.0112706.s006(TIF) Figure S7.adk transcript and gene copies over time under conditions of light/ dark cycling, temperature flux, and continuous darkness.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 derived from triplicate qPCRs from biological replicates.doi:10.1371/journal.pone.0112706.s007(TIF) Figure S8.asl transcript and gene copies over time under conditions of light/ dark cycling, temperature flux, and continuous darkness.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 derived from triplicate qPCRs from biological replicates.doi:10.1371/journal.pone.0112706.s008(TIF) Figure S9.vma3 transcript and gene copies over time under conditions of light/ dark cycling, temperature flux, and continuous darkness.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 derived from triplicate qPCRs from biological replicates.doi:10.1371/journal.pone.0112706.s009(TIF)

Table 1 .
List of candidate reference genes.Thermocycling conditions consisted of 94 ˚C, 3 min; 35 cycles of 94 ˚C, 45 s, 51 ˚C (48 ˚C for act) for 30 s, 72 ˚C, 1 min, and a final extension at 72 ˚C, 7 min.PCR products were visualized by gel electrophoresis to confirm product amplification of the proper size.The resulting amplicons were cloned into the pCR4 vector using the TOPO TA cloning kit (Invitrogen, Carlsbad, CA) and transformed into chemically competent TOP10 E. coli following the manufacturer's instructions.Cells were spread onto Miller's Luria Broth (LB) plates containing 50 mg mL 21 kanamycin (LBkan 50 ) and incubated overnight at 37 ˚C.Colonies were screened for inserts of the expected size via the gene-specific PCRs described above.Colonies producing a band of the expected size were grown overnight in 2 mL LBkan 50 broth.Plasmids were isolated using the QIAprep Spin Miniprep kit (Qiagen, Valencia, CA) and sequenced using the M13 forward and reverse primers by GeneWiz, Inc. (Germantown, MD).
[18]mers from[25].doi:10.1371/journal.pone.0112706.t001bettermodelsfor mRNA quantification due to higher sensitivity, increased quantification range, higher reproducibility, and is more stable than the RNA calibration curve[18].PCRs were conducted in 25 mL 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 mL with nuclease-free water.

Table 3 .
Descriptive statistics and stability values of candidate reference genes under all-dark conditions (n510).

Table 4 .
Descriptive statistics and stability values of candidate reference genes under temperature fluctuations (n512).

Table 5 .
Descriptive statistics and stability values of candidate reference genes under light/dark cycling, all-dark growth, and temperature fluctuations (n530).