Identification of Sexually Dimorphic Genes in the Neonatal Mouse Cortex and Hippocampus

2.1 Mouse body weight, cortex/hippocampus weight, and RNA yield

For the male and female neonatal mice used in the microarray experiment, there were no sex differences in their body weights (p = 0.912), cortex/hippocampus weights (p = 0.084), or RNA yields (p = 0.526) (Table 1). Similarly, the pups used for the RT-qPCR experiment showed no differences in their body weight (p = 0.571), cortex/hippocampal weight (p = 0.423), or RNA yield (p = 0.571) between the sexes (Table 1).

2.2 Microarray screening of sex-biased gene expression in the neonatal mouse cortex/hippocampus

Among 45,281 probes included in the Illumina gene expression microarray (San Diego, CA), 17,270 (38.1%) were detected as being expressed in the neonatal cortex/hippocampus. With the BLAST search of the expressed probe sequences in comparison with the NCBI transcriptome database (RefSeq Release 59), 10,762 (62.3%) matched a single transcript, 3,783 (21.9%) matched multiple transcripts corresponding to a single gene, 774 (4.5%) matched multiple genes, and 1,951 (11.3%) matched nothing in the RefSeq database. A total of 14,545 probes from the first two groups (84.2% of expressed probes) corresponding to 10,489 genes were further tested for sex differences in their expression. Using the false discovery rate (FDR) and 1.2-fold change cutoffs, 95 probes (0.65% of 14,545 expressed probes) were identified as sexually dimorphic; these corresponded to 90 distinct genes (0.86% of 10,489 expressed genes) with more being male-biased (55 genes, 61%) than female-biased (35 genes, 39%) (χ² = 4.44, df = 1, p = 0.035) (Supplemental Table S1).

A summary heatmap for hierarchical clustering analysis of samples and sexually dimorphic candidate genes is shown in Figure 1. First, based on their expression patterns, 17 samples were grouped into two clusters (Fig. 1, the top cluster tree). Nine male samples were clustered together (the right half) before joining the females, and vice versa for females (the left half). Additionally, within each of the clusters, the samples used in the first assay (M7–9 or F7, 8) were connected with the same sex samples from the second (M1–6 or F1–6) before they associated with each other, indicating a lack of significant variation between the two assays. The lack of substantial differences between the two assays was also confirmed by principal components analysis (PCA; Supplementary Figure S1). Similar to clustering analysis, male and female samples were consistently mapped closer to the same sex samples before intermingling with the opposite sex in the three PCA plots generated by using expression data of the 95 sexually dimorphic probes (Panel A), the 494 sex chromosome probes (Panel B), and the 5,000 probes with the lowest expression noise (Panel C). Furthermore, within each plot, the samples from the first assay were intermingled with the same sex samples from the second assay, confirming again no significant differences between these two assays.

Clustering analysis assigned 95 sexually dimorphic gene probes into two clusters, male-biased in the top half and female-biased in the bottom half (Fig. 1, the left cluster tree). Three Y chromosome probes, Eif2s3y, Ddx3y, and Kdm5d (in blue), were more closely clustered together than any of the other male-biased genes, and this cluster was next joined by Cd151. Six female-biased probes, including three from the X chromosome (Eif2s3x, Kdm6a, and Xist; in red) were also tightly clustered and then joined by a cluster of three X/Y pseudoautosomal probes (two for Erdr1 and one for Mid1; in green) and a probe matching 1110017D15Rik on chromosome 4.

Among the 90 sexually dimorphic candidate genes expressed in the neonatal mouse cortex/hippocampus, eleven of them were located on the X and Y chromosomes (3.6% of the 304 sex chromosome genes detected as expressed) (Table 2). Three Y chromosome genes, Ddx3y, Eif2s3y, and Kdm5d, showed the most significant sex difference in expression (all p < 0.001, ranked 1st, 2nd, and 4th most significant). Unlike Eif2s3y and Kdm5d expressed exclusively in male samples, Ddx3y expression was detected in both males and females although it was strongly male-biased. Manually searching the sequence similarity of the Ddx3y microarray probe in the BLAST, we found it matches not only Ddx3y, but also several autosomal transcripts, such as Zdhhc3 (43 out of 50 nucleotides matching) and Orc4 (40 out of 50 nucleotides matching) at a threshold lower than that used in the automated search, suggesting that detection of Ddx3y in females is due to cross-reactivity of the probe with other genes. The RT-qPCR measurement with Ddx3y–specific primers confirmed Ddx3y transcript expressed in males only (Fig. 2A). Xist showed female-biased expression with the largest fold change between the sexes, likely due to its high expression; the four distinct microarray probes matching Xist transcripts were expressed in all female, but not male samples, with 37 to 516 times higher expression in females as compared to males. Aside from Xist, six other probes corresponding to five other X chromosome genes were detected with significant sex differences, including three male-biased (Ebp, Ftsj1, and Psmd10) and two female-biased genes (Eif2s3x and Kdm6a). Kdm6a had two probes both showing female-biased expression with 1.34 (9th most significant) and 1.44 (11th most significant) times higher in females than males, respectively. Eif2s3x ranked 8th in terms of statistical significance, and its expression was 1.41 times higher in females than males.

Three microarray probes matching two sex chromosome genes present in the pseudoautosomal region (PAR) of the X and Y chromosomes were identified as differentially expressed in the male and female mouse cortex/hippocampus: two matching Erdr1 and one matching Mid1. All three probes showed significantly higher expression in females. The two Erdr1 probes ranked 19th and 26th with the fold changes of 1.53 and 1.46, respectively. The Mid1 probe ranked 74th in terms of statistical significance but 18th in terms of fold change, showing expression of 1.56 times higher in females than in males.

Using the Gene Ontology (GO), we examined the functional annotations of these genes and found them involved in a variety of molecular functions (Table 2). For example, Kdm5d and Kdm6a encode enzymes for histone H3 demethylation on K4 and K27, respectively. Eif2s3x and Eif2s3y are involved in translational initiation. Xist involves dosage compensation by silencing one of the two X chromosomes in female cells. Erdr1 negatively regulates cell proliferation, while Psmd10 negatively regulates apoptosis. Ftsj1 showed no GO functional annotation.

2.3 Functional classification of sexually dimorphic candidate genes

Using GO annotation search, we found 26 out of the 90 sexually dimorphic candidate genes with functional linkages to at least one of three biological processes examined (developmental, sex difference, and steroid hormone processes): 20 linked to developmental processes, 6 involved in sex difference processes, and 2 associated with sex hormone processes (Table 3 and Supplemental Table S2). The 20 genes under the category of developmental processes were further explored with their connections to “cell proliferation” (GO: 0008283), “cell death” (GO:0008219), “nervous system development” (GO: 0007399), “cell morphogenesis” (GO:0000902), “cell migration” (GO:0016477), and “synapse organization” (GO:0050808) as well as the descendants of these GO terms (Supplemental Table S2). For example, Klk8 was found to involve “neuron projection morphogenesis” (GO:0048812), a descendant of “cell morphogenesis” (GO:0000902) (Oka et al., 2001). Furthermore, with the other search for “nervous system development” (GO:0007399), Klk8 was also revealed to be involved in “negative regulation of myelination” (GO:0031642) and “negative regulation of axon regeneration” (GO:0048681) (Terayama et al., 2007).

Six candidate genes involved in sex difference processes were similarly analyzed with three GO terms, “sexual reproduction” (GO:0019953), “reproductive system development” (GO:0061458), and “sex differentiation” (GO:0007548) and their descendants (Supplemental Table S2). For example, Prkdc and Dhcr24 were annotated with “germ cell programmed cell death” (GO:0035234) and “male genitalia development” (GO:0030539), which are the descendants of “sexual reproduction” (GO:0019953) and “reproductive system development” (GO:0061458), respectively (Maser et al., 2007; Wechsler et al., 2003). Of the two genes involved in the category of sex hormone processes, Dab2 was linked to "negative regulation of androgen receptor signaling pathway" (GO:0030521) while Ppargc1b was annotated with "intracellular estrogen receptor signaling pathway" (GO:0030520) (Kressler et al., 2002; Zhoul et al., 2005)

2.4 Validation of dimorphic expression of sex chromosome genes in the neonatal cortex/hippocampus

Sixteen out of the 90 sex-biased candidate genes identified by microarray were selected to validate their expression in the neonatal cortex/hippocampus by RT-qPCR, including three on the Y chromosome (Ddx3y, Eif2s3y, and Kdm5d), three on the X chromosome (Xist, Eif2s3x, and Kdm6a), two in the PAR of the X/Y chromosomes (Erdr1 and Mid1), and eight autosomal genes (Cd151, Dab2, Dhcr24, Klk8, Lphn1, Mea1, Meg3, and Prkdc). Eighteen pairs of oligonucleotide primers for these candidate genes and two reference genes were used in RT-qPCR to measure their relative RNA levels (Supplementary Table S3).

Three Y chromosome genes, Ddx3y, Eif2s3y, and Kdm5d, were confirmed to be expressed solely in the male cortex/hippocampus at levels more than 3,000 times higher in males than in females (p <0.001) (Fig. 2A). In contrast, Xist showed expression solely in females with levels over 2,000 times higher in females than in males (p < 0.001) (Fig. 2B, top panel). Unlike Xist, the other two X chromosome genes, Kdm6a and Eif2s3x, were expressed in both sexes, but expression of these genes in females was 2.71 (1.07 ± 0.13 vs. 0.40 ± 0.04, p = 0.001) and 1.28 (1.04 ± 0.07 vs. 0.81 ± 0.05, p = 0.015) times higher than males, respectively (Fig. 2B). Mid1 and Erdr1 are neighboring genes located on the PAR of the X and Y chromosomes as reported by Dal Zotto et al. (1998) and the NCBI Gene database (http://www.ncbi.nlm.nih.gov/gene/); both were confirmed as showing significantly higher gene expression in females, with fold changes of 1.63 (1.05 ± 0.12 vs. 0.65 ± 0.13, p = 0.041) and 1.42 (1.03 ± 0.06 vs. 0.73 ± 0.07, p = 0.002) times greater in females than in males, respectively (Fig. 2C). All of the eight sex chromosome genes tested by RT-qPCR showed significant sex differences consistent with the microarray data.

2.5 Validation of sexually dimorphic expression of autosomal genes by RT-qPCR

Eight autosomal candidate genes listed on Table 3 were selected to validate their sexually dimorphic expression in the neonatal cortex/hippocampus by RT-qPCR. Five genes were confirmed as sexually dimorphic; of those five, four genes were male-biased (Cd151, Dab2, Klk8, and Prkdc) and one was female-biased (Meg3) (Fig. 3). Similar to microarray data (1.34 fold higher in males than in females), Klk8 expression was 1.25 times greater in males (1.33 ±0.10 vs. 1.06 ± 0.08, p = 0.043) (Fig. 3). Relative mRNA levels of Cd151, Dab2, and Prkdc in the neonatal mouse cortex/hippocampus also showed significant male-biased expression with increased fold changes of 1.24 (1.25 ± 0.10 vs. 1.00 ± 0.06, p = 0.041), 1.36 (1.48 ± 0.15 vs. 1.08 ± 0.10, p = 0.038), and 1.26 (1.30 ± 0.11 vs. 1.03 ± 0.07, p = 0.044) in males as compared to females. Meg3 showed a significant female bias with expression 1.36 times higher in the female cortex/hippocampus than in the male (1.10 ± 0.11 vs. 0.80 ± 0.05, p = 0.024) (Fig. 3).

Expression of the remaining three candidate genes, Dhcr24, Lphn1, and Mea1, was not significantly different between the sexes in RT-qPCR. Dhcr24 and Mea1 showed male-bias and Lphn1 displayed female-bias by microarray. Mea1 showed a trend of male-bias (1.18 ± 0.07 vs. 1.03 ± 0.06, p = 0.113) (Fig. 3). Dhcr24 (F vs. M: 1.03 ± 0.07 vs. 1.00 ± 0.06, p = 0.743) and Lphn1 (F vs. M: 1.02 ± 0.06 vs. 1.11 ± 0.08, p = 0.38) showed essentially even expression between the sexes in the neonatal mouse cortex/hippocampus.

4.1 Animals

Adult male and female C57BL/6J mice (Mus musculus) were purchased from the Jackson Laboratory (Bar Harbor, ME) and housed in the California State University, Long Beach (CSULB) Animal Care Facility. Animals were kept on a 12:12 light:dark cycle (lights on at 0600 h). Food (Teklad Mouse/Rat Diet #7012; Harlan Laboratories, Inc., Indianapolis, IN) and water were provided ad libitum. Adult female mice were individually paired with a fertile male, which was removed during late pregnancy. On the day of birth, male and female pups (n= 24–27 per sex) were sacrificed by rapid decapitation. In male mice, a pigment spot was visible on the scrotum, whereas female pups lacked pigmentation in the anogenital region (Wolterink-Donselaar et al., 2009). Their brains were immediately removed from the skulls and coronally blocked between 2.55 mm and 4.95 mm from the most rostral end (Paxinos et al., 2007). The top half of this brain block, including the cerebral cortex and hippocampus, was dissected out, frozen on dry ice, and then stored at −70°C until processed for RNA extraction. Besides examining the presence of pigmentation (males) in their anogenital regions, the tails of these mice were also collected for sex genotyping to confirm their sexes. All experimental procedures were approved by CSULB Institutional Animal Care and Use Committee (IACUC) and performed according to Association for Assessment and Accreditation of Laboratory Animal Care International (AAALAC) guidelines.

4.2 Sex genotyping

Genomic DNA was extracted from the mouse tails using Sigma REDExtract-N-Amp™ Tissue PCR Kit (St. Louis, MO). As described previously (Gatewood et al., 2006; Tsai et al., 2009; Turner et al., 2000), tail DNA was amplified by PCR for the YMT2/B sequence (Ssty1), a member of the Ssty family present on the long arm of Y chromosome, and the myogenin sequence (Myog). PCR products were separated on a 1.5% agarose gel and stained with ethidium bromide. The presence of Ssty1 and Myog was used as the indicator of male genotype and the positive control, respectively.

4.3 RNA Extraction

Total RNA was extracted from the mouse cortex/hippocampus using RNeasy® Lipid Tissue kits (Qiagen Inc., Valencia, CA) according to the manufacturer's protocol. RNA quality and concentration of individual samples was determined using Bio-Rad SmartSpec™ Plus spectrophotometer (Hercules, CA) as the samples were diluted (1:50) in 10 mM Tris (pH = 7.5). Absorbance at 260 nm (A260), 280 nm (A280), and 320 nm (A320) was measured. RNA concentration was calculated with A260, and purity was determined by the A260/A280 ratio (between 1.9 and 2.1).

4.4 Gene expression microarray screening of sex-biased genes

To mitigate the possibility of female embryos being affected by T from neighboring males in utero, the RNA samples used for microarray were from the litters with the prevalent sex; male samples from the litters with the majority of male pups, and female samples from those with the majority of females. Two microarray assays were performed with the pooled cortex/hippocampus RNA samples. In the first assay, RNA of 9 male and 6 female mice (3 littermates of the same sex per sample) was mixed to create three male (M7–9) and two female (F7 and F8) pooled samples. The second assay was performed with six male (M1–6) and six female (F1–6) pooled samples, each of which was composed of RNA from two non-littermates of the same sex (Fig. 1).

To evaluate purity, pooled RNA samples were run on an agarose gel and visualized by ethidium bromide staining, and the 28S/18S rRNA ratio calculated. The RNA samples were amplified, labeled, and hybridized to the Illumina MouseWG-6 Expression BeadChip (San Diego, CA), performed by Expression Analysis® (Durham, NC). Arrays were then scanned by the Illumina BeadArray Reader and processed by BeadStudio software to produce summary data, which showed no difference in the average numbers of genes detected as expressed in the first and the second microarray assays with the p-value thresholds of 0.01 (10168.8 ± 63.0 vs. 10118.0 ± 60.1, p = 0.316) or 0.05 (12427.2 ± 88.2 vs. 12340.8 ± 69.6, p = 0.25). Our microarray data have been deposited in the NIH Gene Expression Omnibus (GEO) database (GSE54283, http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE54283).

4.5 Microarray data analysis

To calculate the fold changes of expression, raw data were first normalized by background subtraction (the median of the negative control probes). Negative values were scrunched into the range 1 to 2 with the same ordering. Normalized data were further processed by quantile normalization (Dunning et al., 2008). Microarray probes were filtered to match current NCBI RefSeq sequences by performing automated BLAST with BioPython (Cock et al., 2009). A probe sequence had to match at least 48 of its 50 nucleotides to be considered effective, and probes matching nothing in the RefSeq database or multiple gene symbols were excluded from further analysis. Filtered probes were tested for expression with significantly higher signal than background (determined by a set of negative control probes) using a one-tailed, one-sample t -test. Probes with significantly higher signals than background (p < 0.01) in at least six of the seventeen total samples were considered expressed. Expressed probes were further tested for a significant difference in expression between the sexes by two-tailed, two-sample t -test on log-transformed expression measurements. Raw p-values were adjusted by false discovery rate (FDR), and an FDR cutoff of 0.05 was used to determine significance (Benjamini and Hochberg 1995). Probes were also filtered for a fold change of at least 1.2 between the sexes. Analysis of the microarray was performed in Python using statistical tools from the SciPy toolkit and storing the data in a MySQL database (Oliphant 2007).

Hierarchical cluster analysis was performed in R with average linkage and a correlation metric to reveal the close relationships among the samples and sexually dimorphic candidate genes. Heatmaps were generated in R using the heatplot function (Eisen et al., 1998). PCA was also performed in Matlab to examine the relationship between the samples based on expression levels of three sets of expressed microarray probes: (1) 95 sexually dimorphic gene probes, (2) 494 X and Y chromosome probes, and (3) 5,000 probes with the smallest mean bead standard errors (noises).

The list of sexually dimorphic candidate genes was searched for the given gene ontology (GO) keywords (Supplementary Table S2) by custom code written in Python to identify annotations of genes with either the listed keywords or any of their descendants in the ontology tree. Annotations were downloaded from Mouse Genome Informatics (MGI) (ftp://ftp.informatics.jax.org/pub/reports/index.html#go) and tree structure was downloaded from the Gene Ontology project (http://geneontology.org/GO.downloads.ontology.shtml).

4.6 cDNA synthesis and PCR primer testing

The cDNA samples (n=18 per sex) were prepared from 1 µg extracted RNA in 20-µl reactions using the Bio-Rad iScript™ cDNA Synthesis Kit. The reverse transcription was performed at 25 C for 5 min, 42 °C for 30 min, and 85 °C for 5 min in a Bio-Rad MyCycler™ (Hercules, CA) thermocycler.

We first looked for PCR primers of selected genes in the literature and in PrimerBank (http://pga.mgh.harvard.edu/primerbank/), and new primers were designed if they were not found. Primers were designed and then checked for suitability using the Primer3 program (http://primer3.wi.mit.edu/), followed by checked for sequence specificity using NCBI’s Primer-BLAST program (http://www.ncbi.nlm.nih.gov/tools/primer-blast/) against RefSeq RNA. All primers found suitable for quantification along with the annealing temperature used and size of the resulting product are listed in Supplementary Table S3.

Primers purchased from Eurofins MWG Operon (Huntsville, AL) were tested in RT-PCR assays using GoTaq® Hot Start Green Master Mix (Promega, Madison, WI) according to the manufacturers’ protocol for 40 cycles. Initial tests used an annealing temperature of 2 °C below the lower melting temperatures of the primers as indicated by Eurofins MWG Operon. PCR products were separated on a 1.5% agarose gel and visualized by ethidium bromide staining. All primers used in this study produce single PCR products with the correct size in the positive control as calculated, but no product in the negative control.

Dilution cDNA stocks were prepared with 8 serial dilutions ranging from 1:8 to 1:1,024. For each qPCR reaction, 4 µl of cDNA, 1 ul each of forward and reverse primers (840 nM), and 6 ul of ABsolute™ QPCR SYBR® Green Mix (Waltham, MA) were mixed and performed on an Agilent Stratagene® MX3000P™ qPCR system with MxPro QPCR software (Santa Clara, CA). qPCR amplification were conducted with an initial 15 min denaturation/polymerase-activation step at 95°C, 40 cycles of 15 s at 95°C, 1 min at the annealing temperature, and 30 s at 72°C, and a dissociation step with 1 min at 95°C, 30 s at 60°C, and 30 s at 95°C. Cycle of threshold (Ct) were analyzed by linear regression with log base 2 of the dilution as the predictor and Ct of the reaction as the response. All primers we used were suitable for quantification with the average R² of the regression of 0.98, Ct values at the 1:16 dilution above 30, and the dissociation curve with single peaks.

4.7 Validation of dimorphic expression of sex-biased candidate genes by RT-qPCR

RT-qPCR for gene expression quantification was conducted under the same conditions as described in the primer dilution tests above, with 4 µl of 1:16 diluted cDNA of the male and female neonatal cortex/hippocampus (n=18 per sex). These samples were randomly selected from all litters; 17 male samples used for RT-qPCR came from male-biased litters, compared to one male from a female-dominated litter. Thirteen female pups used were from female-biased litters, compared to five from male-biased litters. Two replicate reactions for each sample were run and their Ct values averaged before further analysis. For a given sample unnormalized relative expression values for gene expression were calculated using Pfaffl’s method against a female baseline for all cases except Y chromosome genes; for Y chromosome genes, a male baseline was used (Pfaffl 2001). Six reference/housekeeping genes were tested to determine suitability for normalization of relative gene expression. Relative expression values of the target gene were divided by the geometric mean of the relative expression values of the two most stable reference genes (Actb and Rpl13a) to yield normalized relative expression values (Vandesompele et al., 2002). Values more than two standard deviations from their group mean were trimmed to two standard deviations from the group mean in the same direction.

4.8 Statistical Analyses

Aside from the microarray data, the data from individual groups were first tested for normality using Anderson-Darling test, and if passed, the data were subjected to two-tailed, two-sample t -test. In the case of significant deviations from normality, measurements were log transformed, followed by t -test. If transformation still failed to normalize the data, the data were subjected to non-parametric Mann-Whitney test.