Selection and evaluation of reference genes for analysis of mouse (Mus musculus) sex-dimorphic brain development

Sample collection

Whole embryonic brain tissue (overlying epidermis and cranial facial tissues were removed) and tail tips were collected from timed stages (E11.5 to 18.5) pregnant mothers. RNA was extracted from brain tissue using Purelink RNA mini kit (Ambion, USA) according to manufacturer’s instructions. DNA was isolated from tail tips using 0.2 mg/mL (final concentration) Proteinase K (New England Biolabs, Ipswich, MA, USA) and then added to DirectPCR Lysis reagent 102-T (Viagen Biotechnologies, Los Angeles, CA, USA). Tail tips were incubated overnight at 55°C and followed by heat inactivation of the proteinase K at 85°C for 45 min. Samples were centrifuged at 14,000 g for 1 min to pellet cell debris and 2 µL of each sample was used for sexing PCR. The RNA for each time point and sex was a pooled sample, with a minimum of three separate biological samples collected for each condition (from separate litters). Total RNA was quantified with a Nanodrop (ThermoFisher, Waltham, MA, USA) and purity assessed using the 260/230 and 260/280 ratios. A ratio of ∼2 was taken as acceptable for pure RNA.

Sexing of embryos by PCR

To sex the embryos, the Sry gene was amplified using primers listed in Table S1. Following PCR, products were run on a 2% agarose gel, a band appears at approximately 380 bp, indicating the presence of the Y-chromosome for male embryos.

DNase treatment and reverse transcription

To remove genomic DNA present in the sample, 1 µg RNA was added to 1 µL DNase and incubated at 37 °C for 40 min. A phenol/chloroform extraction was carried out followed by ethanol precipitation to purify sample. Reverse transcriptase was performed using iScript (Bio-Rad, Hercules, CA, USA) according to manufacturer’s instructions.

RT-qPCR

Oligonucleotide primers (Table S1) were designed and ordered from Integrated DNA Technologies (IDT, Newark, NJ, USA). Each reaction contained 10 µL SYBR Select Master Mix (ThermoFisher), 1.25 µL 20 pmol forward and reverse primers, 6.75 µL water and 2 µL cDNA. All reactions also included a no reverse transcription control and each reaction was carried out in triplicate. RT-qPCR was carried out in the Stratagene Mx3000p (Agilent Technologies, San Diego, CA, USA) under the following conditions: denaturing at 50°C for 2 min, annealing 96°C for 2 min, followed by 40 cycles of amplification (96°C 15 s, 60°C 15 s and 72°C 1 min). A final cycle for melt curve analysis was included for every qPCR plate (one cycle: 95°C 1 min, 55°C 30 s and 95°C 30 s). The threshold was automatically set by the MXPro program (v. 4.10), and the Ct value for each sample was calculated from the average of three technical replicates.

Data analysis

Data from RT-qPCR was then analyzed using three algorithms to calculate the most stable or best suited reference gene; NormFinder (Andersen, Jensen & Orntoft, 2004), GeNorm with SLqPCR R-based package (Hellemans et al., 2007) Bestkeeper (Pfaffl et al., 2004) and RefFinder (deltaCt method, (Silver et al., 2006; Xie et al., 2012) software). All data was analyzed according to the program instructions.

The geometric mean was calculated for each reference gene by assigning a number (1–10) for top ranked genes for each algorithm. The following equation was used to calculate the geometric mean for each data set: ³√NormFinder*GeNorm*BestKeeper*deltaCt. The reference genes were then ranked again based on their lowest geometric mean.

Analysis of sex specific gene expression was calculated using the following equation: ΔCt = 2^{−(Ctgene-Ctreference)} with the average Ct value for the top two ranked reference genes. Statistical testing used was either the unpaired Student’s t-test when comparing two time-points or a 2-way ANOVA when comparing multiple time-points across time (Tukey’s multiple comparison test).

Selection of candidate reference genes and mRNA transcript levels

To ensure accurate analysis of gene expression, ten candidate reference genes (Table 1) were investigated to determine how stable each gene is for sex and time-point of development. The candidate reference genes selected included those genes commonly used in RT-qPCR experiments (such as Gapdh and Actb) and genes were stably expressed in multiple mouse adult tissues, but had not been tested for suitability as reference genes with mouse embryo tissues (Kouadjo et al., 2007). RT-qPCR was carried out for brain tissue samples across three time points: E11.5, E12.5 and E15.5 using both male and female samples. These time points were chosen for analysis as the sex determination in mouse occurs between E11.0 and E12.0. During this 24 h period, considerable changes also occur with respect to neuronal development, in particular, the formation of the primary brain vesicles (Stiles & Jernigan, 2010). By E12.5, cellular proliferation is increased resulting in the expansion of neuronal precursors and formation of cortical layers (Finlay & Darlington, 1995). Neuronal differentiation, axonal branching and synaptogenesis is taking place around day 15.5 (summarized in Fig. 1) (Sur & Rubenstein, 2005) The raw Ct values for each reference gene across time and each sex are plotted in Fig. 2.

Determination of best-suited reference genes for male and female developing brain

Four algorithms: NormFinder, GeNorm, BestKeeper and comparative deltaCt were employed to determine the best-suited reference genes for either specific developmental stage or between sex over time. Female and male samples were combined to identify the most stable reference gene across time. At each developmental stage (E11.5, E12.5 and E15.5) female and male samples were compared against each reference gene to determine the best reference gene between both sexes.

The raw Ct values for each gene comparing sex-expression at each time point are shown as box and whisker plots (Fig. 2) with minimum and maximum values indicated at each time point. Ct values for male samples ranged between Ct values of 15–25.5 and female samples between 15–28 (Table S3). Across all time points, the reference gene with least amount of variability for male samples is Actb with a mean Ct value of 16.36 (±1.178 SD) while the highest is Eif3f with a mean Ct value of 30.06 (±5.08 SD). In females, Actb mean Ct value is 17.9 (±2.22 SD) and Eef2 has a mean Ct of 22.53, while this Ct value was not the highest, the SD had the highest amount of variability of ±4.54 SD.

NormFinder Excel based add-on is an algorithm created by Andersen, Jensen & Orntoft (2004) that ranks a set of candidate reference genes using stability values according to the variation in expression across samples and between groups. A low stability value represents a reference gene with the most stable expression in a given sample set. NormFinder analysis identified, that across all developmental stages of development (Fig. 3A, Table S4), the reference gene with the lowest stability value was RpL37 (1.211), the best combination of genes were Hprt1 and RpL38 (0.559) in the male embryonic brain. In female samples, Hprt1 (0.929) was the most stable gene across all time points, the best two gene combination was Sdha and RpL38 (stability value of 0.640) (Fig. 3A). NormFinder identified the lowest stability value at E11.5 to be Actb (0.309) (Fig. 3B), however Actb and Sdha (0.250) genes gave the best combined stability value (0.250) for RT-qPCR analysis. Following sex determination of the embryo, RpL38 (0.253) was the most stable reference gene at E12.5, with a combination of RpL38 and RpL37 (0.189) for use of two reference genes (Fig. 3B). RpL37 (0.145) is the most stable reference gene at E15.5 however for analyses with two reference genes then Pgk1 and Eef2 (0.430) were recommended (Table S3).

BestKeeper software determines the ideal reference genes out of a number of candidates based on the standard error (±CP) of Ct values and combines them into an index using repeated pair-wise correlation analysis (Pfaffl et al., 2004). Across all developmental stages, the BestKeeper identified Pgk1 (SE: ±1.064) and Sdha (SE: ±1.004) as the best reference gene in the male sample, while Gapdh (SE: ±1.639) and Sdha (SE: ±1.441) to be ideal for female samples (Fig. 3C). However at specific developmental stages, the ideal reference gene differs. At E11.5 RpL38 (SE: ±0.68) is ranked the top reference gene, while Pgk1 (SE: ±0.522) and Hprt1 (SE: ±1.762) are the most stable reference genes for E12.5 and E15.5 respectively (Fig. 3D, Table S5).

The SLqPCR R-based package uses the GeNorm algorithm (Hellemans et al., 2007) which determines the most stable reference gene from a set of samples by calculating the geometric mean of each reference gene in a stepwise calculation. A low M value indicates stable expression of the gene: for homogenous samples a value below 0.5 indicates an unstable reference gene in the samples analyzed, whereas for heterogeneous tissues the mean M value below 1 is acceptable (Hellemans et al., 2007). The best stability value for the female samples across all time points was Sdha and Gapdh with a stability M value of 0.591 while the top three most stable reference genes in the male samples were Pgk1, and ActB with a mean M stability value of 0.281 (Fig. 4A, Table S6).

At specific developmental time points, Sdha and Pgk1 (M of 0.2470) were calculated as the most stable reference genes for analysis of male and female samples at E11.5. Following sex determination in both male and female samples, Actb and Gapdh (M of 0.346) gene are the most stable (Fig. 4B). During increased neuronal proliferation at E15.5, Actb and Pgk1 (M of 0.719) gave the lowest M values, indicating these are reliable reference genes according to GeNorm analysis given at this stage the developing brain is heterogeneous (Table S6).

A forth method, termed the comparative deltaCt method (Silver et al., 2006) was also used to compare reference genes. This uses a method similar to that of GeNorm, determining the variation of gene expression between paired putative reference genes (the deltaCt) within each sample. Variable deltaCt values for each pair between multiple samples, resulting in a high standard deviation when the average deltaCt is calculated, means either gene or both are not stably expressed. Pairs at E11.5 were Pgk1 and Sdha (1.98/1.87), Actb and Pgk1 (1.98/1.82). At E15.5 there was much more variation between reference genes, best pair were Actb and Sdha (4.7) (Fig. 4D; Table S7) When considering brain development overtime, for the male samples Pgk1 and Sdha (2.1) were the best pair, whereas for female samples, Hprt1 and Sdha (2.7) had the lowest stand deviation of deltaCt values (Fig. 4C).

Geometric mean of each reference gene

Each reference gene was ranked from most stable to least stable, as determined by the four stability calculations described above, and assigned a number from 1 (most stable) to 10 (least stable) (Table S8). The geometric mean was taken for each set of rankings to calculate an overall ‘best’ reference gene for each condition (Tables 2 and 3). The recommended reference genes across all stages in both male and female samples are Sdha and RpL37 (GeoMean: 1.18 and 2.94 respectively). Across all stages of embryonic development in males, the recommended reference genes are Actb and Sdha (GeoMean: 2.73 and 1.96) while in the female brain, the most stable set of reference genes are Sdha and Gapdh (GeoMean 1.2 and 2.34).

For RT-qPCR analyses at specific developmental time points regardless of sex (Table 3), Sdha and Pgk1 (GeoMean: 1.56 and 1.86) are best suited for analysis at E11.5. Following determination of sex in the embryo at E12.5, RpL38 and Sdha with GeoMean values of 1.68 and 1.73 are recommended. Finally, RpL37 and Actb (GeoMean 1.86 and 1.86) are the best suited for analysis of E15.5 brain tissue.

In comparison to a careful sex or stage approach above, we also pooled all the data together (both sexes and all stages). Overall, this approached produced stability values that were high (indicating a great variation between samples) for GeNorm and BestKeeper packages (Tables S5 and S6). The geometric mean of ranked references genes for a pooled sample approached suggested that Sdha and RpL37 (Tables 2 and 3), despite RpL37 often being ranked a poor choice if we were just consider the stabilities values for E11.5 and E12.5 individually (Table S8 and S3).

Confirmation of sex-specific gene expression

Previous micro-array data showed sex-specific expression of a number of genes within the embryonic brain (forebrain) at E10.5 (Dewing et al., 2003). Wingless-type MMTV integration site family, member 10B (Wnt10b) was shown to be up-regulated in male embryos, with 1.7 fold increase in expression at E10.5 with respect to the female brain. Cytochrome P450,7b1 (CYP7B1) was expressed 2.1 fold higher in the male brain, whereas expression of X-inactive specific-transcript (Xist) was confirmed significantly higher in the female embryo with a fold change of 18.5 compared to male samples. However, these differences in gene expression were not studied over further time points, when neurogenesis has commenced. To extend the findings from Dewing et al. (2003), regarding sex-dimorphic expression of these genes, the top two ranked reference genes at each time point from this study are used to compare against a selection of three genes in male and female samples across embryonic brain development.

To normalize the expression of our three genes of interest, Sdha and Pgk1 oligonucleotide primers were used at E11.5, RpL38 and Sdha were used for analysis for E12.5. For analysis of E15.5, the reference genes Actb and RpL37 were used. Sex-dimorphic expression of Xist was detected at all three time points with increased expression in female brain compared to the males (Fig. 5A). CYP7B1 and Wnt10b also showed significant sex-dimorphic expression at E11.5 (P < 0.05), with higher expression observed in the male brain at this stage (Figs. 5C and 5E; 3.5-fold and 5-fold respectively). Additionally, following sex determination, Wnt10b and CYP7B1 expression was not significantly different between the sexes using these stage-specific reference genes (Figs.5C and 5E). In contrast, when pooling all data from sex and age together the top two ranked genes were RpL37 and Sdha. When using these two genes to normalize expression data there is a large variation in the expression of the three genes, especially at E12.5, between biological replicates (Figs. 5B, 5D and 5F). No significant differences were observed in gene expression between male and female samples, when using Pgk1 and Sdha as reference genes with the exception of Xist expression at E12.5 (Figs. 5B, 5D and 4F).

When studying gene expression changes over time, two slightly different sets of reference genes were ranked highest, depending upon the genotype of the sample (Table 2). We normalized gene expression at each time point for female samples (with Sdha and Gapdh) and for male samples (with Pgk1 and Sdha) to study how expression of these genes changes with respect to developmental age (Fig. 6). In comparison, when RpL37 (one of the poorer stable genes for female and male time points) and Sdha where used as reference genes for normalization of the all data points (top two genes when reference gene data was pooled Table 2), this introduced so much variation between samples, any analysis would be inconclusive (Fig. S1).

CYP7B1 mRNA expression declines between E11.5 and E12.5 in both sexes but its expression significantly increases again in the female only by E15.5 (Fig. 6A, P < 0.05). Xist transcript expression in the female brain also increased between E11.5 and E15.5 (following sex-determination window (Fig. 1)) (Fig. 6B; P < 0.001). Gene expression of the Wnt10b gene in the developing brain was low at all time-points and did not change significantly overtime in the female staged embryo samples (Fig. 6C). However, the higher expression observed in the male embryonic samples decreased also between E11.5 and E12.5, to similar expression levels found in XX samples (Fig. 6C).