Skip to main content

Advertisement

Sexual dimorphism in brain transcriptomes of Amami spiny rats (Tokudaia osimensis): a rodent species where males lack the Y chromosome

Article metrics

Abstract

Background

Brain sexual differentiation is sculpted by precise coordination of steroid hormones during development. Programming of several brain regions in males depends upon aromatase conversion of testosterone to estrogen. However, it is not clear the direct contribution that Y chromosome associated genes, especially sex-determining region Y (Sry), might exert on brain sexual differentiation in therian mammals. Two species of spiny rats: Amami spiny rat (Tokudaia osimensis) and Tokunoshima spiny rat (T. tokunoshimensis) lack a Y chromosome/Sry, and these individuals possess an XO chromosome system in both sexes. Both Tokudaia species are highly endangered. To assess the neural transcriptome profile in male and female Amami spiny rats, RNA was isolated from brain samples of adult male and female spiny rats that had died accidentally and used for RNAseq analyses.

Results

RNAseq analyses confirmed that several genes and individual transcripts were differentially expressed between males and females. In males, seminal vesicle secretory protein 5 (Svs5) and cytochrome P450 1B1 (Cyp1b1) genes were significantly elevated compared to females, whereas serine (or cysteine) peptidase inhibitor, clade A, member 3 N (Serpina3n) was upregulated in females. Many individual transcripts elevated in males included those encoding for zinc finger proteins, e.g. zinc finger protein X-linked (Zfx).

Conclusions

This method successfully identified several genes and transcripts that showed expression differences in the brain of adult male and female Amami spiny rat. The functional significance of these findings, especially differential expression of transcripts encoding zinc finger proteins, in this unusual rodent species remains to be determined.

Background

The undifferentiated gonad and brain are programmed to be either male or female in most sexually reproducing species. Brain sexual differentiation in therian mammals is generally influenced by a surge in gonadal steroid hormones during embryo development followed by a second surge later in adulthood. This paradigm is termed “organization-activational programming” and occurs in a sex-specific manner [1,2,3]. The organizational period is typified by spikes in androgens and/or estrogens during development and results in permanent masculinization and defeminization of the brain in males [4,5,6,7]. Masculinization of several brain regions, including the hippocampus, is modulated by aromatization of testosterone to estrogen [4,5,6,7]. Manifestation of sex-dependent behaviors, especially male-specific traits, requires a second spike in steroid hormones at adulthood (“activational period”) [8,9,10,11,12,13]. The salient question remains as to the potential contribution of sex chromosome-associated genes in guiding brain sexual differentiation in therian mammals.

Sex-determining region Y (Sry) resides on the Y chromosome in therian mammals, which includes placental mammals and marsupials [14,15,16,17,18]. This gene serves as a transcriptional activator modulating testis development, and thereby, it was coined the testis determining factor (TDF) [19,20,21]. SRY includes an HMG box as a DNA-binding domain [22] that activates in the testes Sry-related box 9 (Sox9) [23], cerebellin 4 precursor gene (Cbln4) [24], transcription factor 21 (Tcf21, former also called Pod1) [25], and neurotrophin 3 (Nt3) [26]. However, the role of Sry in brain sexual differentiation is uncertain.

The four core genotype (FCG) mouse model has been used in elucidating the neural effects of Sry relative to other Y chromosome associated genes [27,28,29,30,31,32,33,34]. In this model, Sry is deleted from the Y chromosome and re-inserted as a transgene on an autosomal chromosome in both XX and XY chromosome bearing mice [27]. The resulting breeding scheme gives rise to the FCG mouse model as offspring from these mated pairs can be one of four different genotypes: XXSry(karotypically and gonadally female, naturally lacking Sry), XXSry+ (karoytypically female but gonadally male due to presence of autosomal Sry transgene), XYSry (karyoptically male but gonadally female due to deletion of endogenous Sry gene), and XYSry+ (karyotypically and gonadally male) [28]. Results to date with this model reveal that Sry and other sex chromosome associated genes interact with steroid hormones to affect neurobehavioral programming. Gonadectomized XX females consume more food and show increased adiposity relative to XY mice [29]. On the other hand, intact XXSry and XYSry mice have increased activity levels, consume less food, and show enhanced anxiety-like behaviors [30]. Social and parenting behaviors are also influenced by Sry and other Y chromosome associated genes [32, 35]. This model indicates Sry regulates neural expression of progesterone receptor (PR) in the anteroventral periventricular nucleus (AVPV), medial preoptic area (MPOA), and ventromedial nucleus, gamma-aminobutyric acid (GABA)/serotonin/dopamine-related genes within the frontal cortex, and growth hormone (Gh) expression in the hypothalamus [27, 33, 34].

To ascertain the potential role of the Y chromosome and Sry in regulating brain sexual differentiation, a better approach though would be to examine the brain transcriptome profile in males and females of therian mammals where the males lack a Y chromosome/Sry and both sexes possess an XO system. Such is the case with two Tokudaia (Amami [Ryukyu] spiny rat- Tokudaia osimensis and Tokunoshima spiny rat- T. tokunoshimensis) and two Ellobius (Transcaucasian mole vole- E. lutescens and Zaisan mole vole- E. tancrei) species [36,37,38,39,40,41]. Both Tokudaia species are critically endangered, and thus, it is vital to begin to understand how gonad and brain sexual differentiation occurs in them. In Amami spiny rats, two possible genes that might guide testes formation are chromobox protein homolog 2 (Cbx2), which acts upstream of Sry, and Er71, also called Ets variant 2 (Etv2), which is upregulated by Sry in species containing this gene [42, 43]. However, it is not clear which genes or transcript isoforms show sex differences in this Tokudaia species. Based on the current status of Amami spiny rats, it is not permissible to obtain brain samples from embryonic or neonatal individuals. Thus, it is not possible to determine which gene(s)/transcript(s) might guide gonad and potentially brain sexual differentiation in these species. Sex differences in brain expression remain evident in adult mice, birds, and humans [44,45,46,47,48]. To begin to understand, neural sex differences in Amami spiny rat, the current studies sought to determine which genes/transcript isoforms show sexually dimorphic patterns of expression in the brain of this species. An ancillary goal of this study was to develop the first reference transcriptome library for Amami spiny rat that may aide in future studies designed to identify putative sex-determining gene(s)/transcript(s) in this species. RNAseq analyses was performed on whole brain samples from male and female adult individuals that died accidentally as part of a government-approved mongoose eradication project on the island in which this species inhabits. The overarching hypothesis was that there would be considerable sexually dimorphic differences in the brain transcriptomic profile with signature patterns identified in each of the sexes.

Results

General characterizations

Table 1 lists the summary statistics of the transcriptome assembly for the Amami spiny rat brain samples. A total of 569,369 transcripts were assembled with a mean length of 839.88 nucleotides (nt). 116,192 transcripts were at least 1000 nt in length. We used a BUSCO analysis to gauge quantitatively the completeness of the assembly as compared to a database of single copy vertebrate orthologs sequences [49]. The Amami spiny rat transcriptome assembly contained 76.4% of the 4104 single copy vertebrate BUSCO orthologs. 29.4% of the orthologs were also single copy in the Amami spiny rat assembly, whereas 47% were duplicated. 16.1% were fragmented and 7.5% were missing. These numbers are similar to other assembled transcriptomes for samples collected from non-inbred, wild vertebrate populations [50,51,52,53]. Full details on the transcriptome library generated for this species are included on https://www.ncbi.nlm.nih.gov/bioproject/474959. Relationships between samples were visualized by using principal component analyses (PCA). Upon plotting the PCA results, three samples were obvious outliers: male sample #1, female sample #2 and female sample #3, as evidenced by the fact these samples were at a significantly different location from other points in the PCA plot (Fig. 1). These samples were removed from further expression analysis.

Table 1 Transcript Assembly Metrics for the Transcriptomic Library Generated from the brain of male and female Amami spiny rats
Fig. 1
figure1

PCA analysis. Three samples were obvious outliers: male sample #1, female sample #2 and female sample #3, and thus, these samples were removed from further RNAseq expression analysis

Differential expression of genes in the brain of male and female Amami spiny rats

The RNA-Seq data was aligned to the sequences in the transcriptome assembly, described above, using RSEM (48), which also estimates the expression levels of transcripts. We further processed gene expression levels using tximport [54]. Differential expression analysis used DESeq2 [55]. Only transcripts having a sum of counts from all samples greater than 25 were analyzed. The genes that displayed differential expression (q ≤ 0.05) in both analyses were examined further (Additional file 1). Volcano plot analyses revealed that several genes showed differential expression in the brain of male compared to female Amami spiny rats (Fig. 2). As shown in Tables 2, 34 genes were upregulated in the brain of male compared to female Amami spiny rats. These include such genes as lumican precursor (Lum), seminal vesicle secretory protein 5 precursor (Svs5), cytochrome P450 B1 (Cyp1b1), and a few zinc finger protein-like genes. In contrast, only four genes were downregulated in males compared to females, including serine (or cysteine) protease inhibitor, clade A, member 3 N isoform X1 (Serpina3n), V(D)J recombination-activating protein 2 (Rag2), antigen KI-67 isoform X3 (Mki67), and B2 bradykinin receptor isoform XI (Bdkrb2) (Table 3).

Fig. 2
figure2

Volcano plot analyses of differentially expressed genes. The volcano plot represents the relationship of each gene’s log2 fold change vs -log10 pvalue. Orange points represent expression ratios greater than 2X (either up or down), red points represent genes whose FDR-corrected p-values are < 0.05 and green points are genes where both properties occur (fold-change > +/− 2X and FDR < 0.05)

Table 2 Top genes upregulated in the brain of male compared to female Amami spiny rats
Table 3 Top genes downregulated in the brain of male compared to female Amami spiny rats

Differential expression of transcripts in the brain of male and female Amami spiny rats

Next, we considered brain transcripts that may show differential expression in the brain of male vs. female Amami spiny rats and vice versa. With this approach, 432 transcripts were upregulated in males and 508 transcripts are downregulated in males (volcano plot shown in Fig. 3; Additional file 2). A subset of transcripts upregulated in males is listed in Table 4. As shown, several transcripts from epigenetic-associated genes were upregulated in males, such as Sin3-Associated Protein P30-Like (Sap30l), SET domain and mariner transposase fusion gene (Setmar), methyl-CpG-binding domain protein 4 (Mbd4), and histone deacetylase 2 (Hdac2). Transcripts from genes encoding several zinc finger proteins were also increased in the brain of males compared to females. Some examples of these genes include zinc finger protein X-linked (Zfx), zinc finger and BTB domain-containing protein 44 (Zbtb44), zinc finger protein 708 (Zfp708), LOC102640673, zinc finger protein (Zpr1), and zinc finger protein 655 (Zfp655). Transcripts from genes associated with male reproduction were also elevated in the brain of Amami spiny rat males, such as testis-expressed protein 9 (Tex9) and spermatogenesis Associated 7 (Spata7).

Fig. 3
figure3

Volcano plot analyses of differentially expressed transcripts. The volcano plot represents the relationship of each transcript’s log2 fold change vs -log10 pvalue. Orange points represent expression ratios greater than 2X (either up or down), red points represent transcripts whose FDR-corrected p-values are < 0.05 and green points are transcript where both properties occur (fold-change > +/− 2X and FDR < 0.05)

Table 4 Select transcripts upregulated in the brain of male compared to female Amami spiny rats

Examples of transcripts that were downregulated in males (or upregulated in females) are listed in Table 5. Unlike those that were upregulated in males, these transcripts are derived from genes representing diverse functions, although select examples of epigenetic-associated genes are included on this list, e.g. lysine acetyltransferase 6A-like (Kat6al), and zinc finger proteins (e.g., zinc finger protein 410 [Zfp410]). The list of the top 35 transcripts upregulated in females includes the associated genes for, mitofusin 1 (Mfn1) and mitofusin 2 (Mfn2), which encode for GTPases in the outer membrane of the mitochondria. The list also includes transcripts from the genes cullin 1 (Cul1) and cullin-3 (Cul3), which encode hydrophobic proteins that provide a scaffold for post-translational modification of cellular proteins by ubiquitin ligases (E3). Transcripts derived from several ubiquitin ligase genes were upregulated in the brain of female Amami spiny rats (ring finger protein 220 [Rnf220], SNF2 histone linker PHD RING helicase [Shprh], Praja ring finger ubiquitin ligase 2 [Pja2], Siah E3 Ubiquitin Protein Ligase 2 [Siah2], Mahogunin ring finger 1 [Mgrn1], and ring finger protein 25 [Rnf25]).

Table 5 Select transcripts down-regulated in the brain of male compared to female Amami spiny rats

Gene interaction metanalyses and pathway analyses

We used STRING [56] to identify potential interactions between differentially expressed genes. Based on the genes identified to be upregulated in males, two main STRING clusters were identified with one cluster comprised of genes such as myosin light chain 1 (Myl1), myoglobin (Mb), triadin (Trdn), myosin regulatory light chain 2 (Myl2), cysteine and glycine rich protein 3 (Csrp3), and small muscle protein X-linked (Smpx) in the center (Fig. 4). The second cluster includes fibromodulin (Fmod), osteomodulin (Omd), osteoglycin (Ogn), decorin (Dcn), collagen type III alpha 1 chain (Col3a1), and Lum. Pathways that were enriched for genes upregulated in males include collagenesis, myofibril formation, and proteinaceous extracellular matrix (Additional file 3). The four genes (Mki67, Bdkrb2, Rag2, and Serpina3n) upregulated in females do not share any similar pathways, and thus, do not form a distinct cluster.

Fig. 4
figure4

STRING pathway analyses diagram for genes upregulated in males. Based on genes identified to be upregulated in males, two main clusters are identified with one cluster comprised of genes such as Myl1, Mb, Trdn, Myl2, Csrp3, and Smpx in the center. The second cluster includes Fmod, Omd, Ogn, Dcn, Col3a1, and Lum. Network nodes represent proteins. Splice isoforms or post-translational modifications are collapsed, i.e. each node represents all the proteins produced by a single, protein-coding gene locus. Node colors are arbitrary. Edges represent protein-protein associations. Associations are meant to be specific and meaningful, i.e. proteins jointly contribute to a shared function; this does not necessarily mean they are physically binding each other. Edge confidence is the relative amount of supporting evidence for the connection between two proteins in the diagram. Thicker lines have more support; thinner lines have less support

When individual transcripts are considered, the primary group that is upregulated in the brain of male Amami spiny rat are zinc finger proteins (Fig.5, yellow highlighted area). Several other smaller clusters are present. Pathways that are enriched in transcripts upregulated in males include Krueppel-associated box and zinc finger pathways (zinc finger zinc finger, C2H2-type/integrase DNA, zinc finger, C2H2, and zinc finger C2H2-like (Additional file 4). Processes associated with upregulated male transcripts include, metabolism, tissue development, regulation of ATPase activity, regulation of several aspects of organelle function, nucleoplasm, microtubule cytoskeleton, and GABA-A receptor complex (Additional file 4).

Fig. 5
figure5

STRING pathway analyses diagram for transcripts upregulated in males. The primary group of transcripts upregulated in the brain of male Amami spiny rat are transcripts encoding zinc finger proteins (yellow highlighted area). Several other smaller clusters are present. Network nodes represent proteins. Splice isoforms or post-translational modifications are collapsed, i.e. each node represents all the proteins produced by a single, protein-coding gene locus. Node colors are arbitrary. Edges represent protein-protein associations. Associations are meant to be specific and meaningful, i.e. proteins jointly contribute to a shared function; this does not necessarily mean they are physically binding each other. Edge confidence is the relative amount of supporting evidence for the connection between two proteins in the diagram. Thicker lines have more support; thinner lines have less support

The transcripts upregulated in females represent diverse pathways (Fig. 6, Additional file 5), including nucleotide, purine ribonucleoside, ribonucleotide, carbohydrate derivative, heterocyclic, zinc ion, ATP, transitional metal ion, cytoskeletal protein, and GABA receptor binding. Additional pathways include hydrolase, nucleoside triphosphate, ATPase, L-amino acid transmembrane transporter, and palmitoyl-CoA 9-desaturase activity.

Fig. 6
figure6

STRING pathway analyses diagram for transcripts upregulated in females. Based on the elevated transcripts, several pathways are affected in females. Examples include nucleotide, purine ribonucleoside, ribonucleotide, carbohydrate derivative, heterocyclic, zinc ion, ATP, transitional metal ion, cytoskeletal protein, GABA receptor binding pathways, hydrolase, nucleoside triphosphate, ATPase, L-amino acid transmembrane transporter, and palmitoyl-CoA 9-desaturase activity. Network nodes represent proteins. Splice isoforms or post-translational modifications are collapsed, i.e. each node represents all the proteins produced by a single, protein-coding gene locus. Node colors are arbitrary. Edges represent protein-protein associations. Associations are meant to be specific and meaningful, i.e. proteins jointly contribute to a shared function; this does not necessarily mean they are physically binding each other. Edge confidence is the relative amount of supporting evidence for the connection between two proteins in the diagram. Thicker lines have more support; thinner lines have less support

qPCR analyses

Due to limited remaining tissue, only three candidate genes shown to be altered in the RNAseq analyses were validated with qPCR. These include two that were upregulated in males: Svs5 and Cyp1b1 and one that was upregulated in females: Serpina3n. The three genes chosen all demonstrated significant adjusted p values and fold change and have previously ascribed roles in guiding sexual differentiation and/or neural function. As shown in Fig. 7, Svs5 expression was considerably upregulated in the brain of males compared to females (p = 0.0005); whereas Serpina3n was significantly down-regulated in males vs. females (p = 0.05), which is identical to the profile detected with RNAseq. While Cyp1b1 showed a tendency to be increased in males, as seen with RNAseq, the qPCR results were not significant.

Fig. 7
figure7

qPCR analyses for select genes in the brain of female and male Amami spiny rats. A) Svs expression. B) Cyp1b1. C) Serpina3n. *p = 0.0005. **p = 0.05

Discussion

When SRY was identified, it was presumed that it regulated sexual differentiation in all eutherian mammals [14,15,16,17,18]. However, those species that defy our expectations reveal the complexity and elegance of sexual differentiation that like other biological processes has been subject to evolutionary forces million years in the making. Such is the case for select species within Tokudaia and Ellobius. Intriguingly, 46, XX human male brothers lacking SRY possess testes but are azoospermic [57]. Overexpression of SOX9 in these SRY-negative men likely stimulates male sexual differentiation [58, 59]. Thirty-eight intersex, XX SRY-negative cases have also been reported in pigs [60].

How sexual differentiation occurs in rodent species where males lack a Y chromosome and SRY remains enigmatic. In Amami spiny rats, two possible candidates have been postulated to regulate testes formation: Cbx2, which acts upstream of Sry, and Er71 [42, 43]. The wave of sexual differentiation affects both the gonad and brain [1,2,3]. While in the latter steroid hormones, testosterone and estrogen, are the guiding factors, SRY might also be involved [61]. In those species that possess SRY, it continues to be expressed in the adult male brain [62, 63]. At the current time, we cannot obtain embryonic or neonatal brain samples from Amami spiny rat to identify putative sex determining gene(s)/transcript(s) in this species. Our prediction at the outset was that similar to mice, birds, and humans [44,45,46,47,48], sex differences would be identified in gene/transcript expression in adult Amami spiny rat, albeit many of these many not necessarily be sex-determining genes. The additional objective of the current work was to establish the first reference transcriptomic library for Amami spiny rat that would likely facilitate future work seeking to identify those gene(s)/transcript(s) compensating for the absence of the Y chromosome and SRY in this species. With these limitations and potential strengths in mind, we analyzed the transcriptomes of the adult male and female Amami spiny rat brains to determine which genes and transcripts show sexually dimorphic patterns of expression.

While certain genes differed in expression between male and female Amami spiny rat, as detailed further below, more prominent changes were discovered when comparing sex differences of individual transcripts. The aggregate characteristics of the differentially expressed transcripts are compelling. While the proportional contribution of each splice variant to the total mRNA for a given gene might be low, such variants might demonstrate different potencies in binding and activating various cellular pathways. Studies with mice and other taxa demonstrate that expression of such transcript variants can result in dramatic sex differences and potentially underpin sexual differentiation [64,65,66,67]. In Amami spiny rat, additional phenotypic studies are needed to validate the significance of these findings and whether they broadly differ across various sexually dimorphic organs, including the gonad, and discrete brain regions.

Several transcripts upregulated in the whole brain of males compared to females encode zinc finger proteins. Additionally, gene set enrichment analyses of the genes from which the upregulated transcripts derive identified several zinc-finger related pathways, including Krueppel-associated box (KRAB). KRAB only (KRAB-O) serves as an SRY-interacting protein [68, 69]. Zinc finger proteins are polypeptides exhibiting sequence-specific, nucleic acid-binding role, and serve as trans-acting molecules they govern cellular growth and differentiation. Zinc finger protein Y-linked (Zfy) is another gene located on the Y chromosome, and while not the TDF, it is essential for normal spermatogenesis [70], although an isolated human case has been reported of normal sexual differentiation in the absence of ZFY [71]. Similar to SRY, this gene continues to be expressed in the brain of adult men [62]. In the Amami spiny rat, Zfy, along with eukaryotic translation initiation factor 2 subunit 3, Y-linked (Eif2s3y) and lysine Demethylase 5D (Kdm5d), have been translocated to the X chromosome [72]. The current brain transcriptome studies detected Zfy, but the gene and individual transcripts showed similar expression pattern in males and females.

A corollary gene is present on the X chromosome, Zfx with male mice harboring a mutation of this gene undergoing normal gonadal and male genitalia sexual differentiation but showing markedly reduced sperm counts [73]. Mutant Zfx female mice possess normal ovaries and external genitalia but show diminished oocyte populations and reduced fecundity culminating in abbreviated reproductive potential. Primordial germ cell populations prior to genital ridge migration are reduced in both XX and XY mice offspring and growth impairments are evident in both sexes. The collective data suggest that Zfx might thus be important in both male and female sexual differentiation in eutherian mammals. In Drosophila, the hermaphrodite (her) gene encodes for a C2H2-type zinc finger protein that demonstrates similar properties as Zfy and Zfx and presumably regulates terminal function in sexual differentiation [74]. In the current studies, one form Zfx was shown to be upregulated in the brain of males. However, several other genes encoding zinc finger proteins were also upregulated in males. Thus, we postulate that upregulation of Zfx acting possibly in concert with other C2H2-type zinc finger transcript forms in male Amami spiny rat might compensate for the absence of SRY and serve as the collective factors stimulating male sexual differentiation.

In support of this hypothesis, other zinc finger proteins exert essential sexual differentiation roles in lower animal species and mice after initial sexual differentiation. For instance, the flexuosa (fle1), encoding aC2H2 zinc finger protein, is an important repressor of female sexual differentiation in Podospora anserina [75]. The byr3 zinc finger protein (byr3) gene encoding a protein with seven finger domains affects sexual differentiation pathways in Schizosaccharomyces pombe [76]. In the Japanese wrinkled frog (Rana rugosa), an increase in testosterone causes sex reversal from female to male and accompanying decreases in zinc finger proteins, zinc finger protein 64 (Zfp64) and zinc finger protein 112 (Zfp112) [77]. In sex-reversed and normal mice, zinc finger protein 35 (Zfp35) is essential for normal spermatogenesis [78, 79]. Zinc finger proteins are elevated in human testicular cell lines overexpressing SRY, further highlighting the potential role of such proteins in maintaining sexual differentiation [80].

Sex differences in transcript isoforms encoding zinc finger proteins and other transcripts might originate due to epigenetic variation. In this aspect, intragenic DNA methylation and histone protein modifications can affect RNA polymerase II activity and thereby stimulate creation of alternative splice forms [81,82,83,84,85]. Heterochromatin protein 1 (HP1) might also play a role in DNA methylation effects on mRNA splicing [86]. DNA methylation and histone protein changes are evident in turtle species demonstrating temperature sex determination [87,88,89,90,91]. In the case of doublesex and mab-3 related transcription factor 1 (Dmrt1), DNA methylation status tightly associates with temperature and may thus be the master male sex determining gene in red-eared slider turtles (Trachemys scripta) [87].

RNAseq analyses revealed other genes and transcripts that were differentially expressed in the brain of males and females. It is not clear if such genes contribute to the initial brain sexual differentiation process or maintaining sex differences. Due to limited tissue, we could only validate three differentially expressed genes that were chosen based on there highly significant adjusted p value and overall fold change: Svs5, Serpina3n, and Cyp1b1. Of these, Svs5 and Serpina3n expression patterns were identical to that observed with RNAseq with an increase in males and females, respectively. Further work is needed to determine the significance of these and other gene expression changes.

The current studies with adult Amami spiny rat yielded less genes that show sex-dependent brain expression than past studies with mice, birds, and humans [44,45,46,47,48]. This difference could be attributed to the fact that from an algorithmic standpoint, potentially high level of intra-sample variability would decrease the number of inter-sample genes whose differential expression fell below a FDR of 0.05. The current work presumably identified those genes that show the greatest differential expression between sample types and are potentially the most important.

One of the primary limitations of this study is that in this endangered species, we could only screen random brain tissues in limited number of adult male and female Amami spiny rats. It is not clear if genes and transcripts guiding initial brain sexual differentiation continue to be expressed. The fact that ZFY and SRY are expressed in the brain of adult men shores up this possibility [62, 63]. Usage of induced pluripotent stem cells (iPSC) generated from this species might be useful in validating the current findings [92, 93]. Future work should compare current RNAseq results in the brain of Amami spiny rat to those obtained in Tokudaia species that possess sex chromosomes, e.g., Okinawa spiny rat (Tokudaia muenninki) [94,95,96], as this approach may reveal those genes/transcripts essential in sex determination in Amami spiny rat. Should the population of these Tokudaia species recover and/or permission be granted in the future to attempt to breed them in captivity, follow-up transcriptome studies should be conducted with discrete brain regions identified to show sex differences in other species, such as the hypothalamus, and with individuals spanning from the embryonic to adult period. As detailed in the methods, one of the other limitations of this study is that it depends upon whole brain samples from animals who died of accidental causes during a mongoose eradication project. In these cases, brain and other tissues were collected several hours to half a day post-mortem. Thus, these studies should be considered a first pass attempt to obtain a reference transcriptome for this species and examine for sex-dependent differences in gene and transcript expression in the brain of adult individuals. A better understanding of underpinning biological mechanisms in this species may aide in recovery efforts to the point where sufficient animals can be bred in captivity and more precise studies in the brain, gonad, and other sex-dependent organs performed.

Conclusions

In conclusion, transcriptomic analyses in the brain of adult male and female Amami spiny rats shows that sex differences are observed in select genes with more being upregulated in males than females. Comparison of individual transcript forms between males and females yielded more differences and potentially might provide mechanistic insight into how male sexual differentiation occurs in this Y chromosome deficient and SRY-negative rodent species. The most attractive candidates are upregulation of a variety of transcript forms encoding zinc finger proteins, including Zfx. Such sex differences in these and other alternative splice forms might originate due to epigenetic modifications. It remains to be determined if any of these identified gene(s)/transcript(s) modulate initial sexual differentiation in this species.

Methods

Animal tissue collection and RNA isolation

T. osimensis is endangered (The IUCN Red List of Threatened Species; https://www.iucnredlist.org/). This species has been protected by the Japanese government as a natural monument since 1972 and Nationally Endangered Species of Wild Fauna and Flora since 2016. With permission from the Agency for Cultural Affairs and the Ministry of the Environment in Japan, tissues were harvested from animals who died accidentally as part of a Japanese government approved mongoose eradication project on the island this species inhabits. A copy of the government approval, which has also been translated into English, to harvest brain and other tissues from animals that die accidentally has been provided to the journal.

The brain and other tissues were harvested and frozen several hours to half a day post-mortem. Brain tissues were scraped by a needle or tweezer through the posterior cranial fossa in order to preserve the skull, which is also rare and needed for taxonomic studies. Several portions of brain tissue were then mixed and frozen together. Frozen male (n = 3) and female (n = 4) brain samples stored at − 80 °C were placed in RNAlater™-ICE Frozen Tissue Transition Solution (ThermoFisher Scientific, Waltham, MA) and stored overnight at − 20 °C per the manufacture’s protocol. The samples were then shipped to Dr. Rosenfeld’s laboratory at the University of Missouri. RNA was isolated from each brain sample with TRIzol Plus RNA Purification Kit (ThermoFisher Scientific). The quantity and quality of the RNA was determined with a Nanodrop ND1000 spectrophotometer (Nanodrop Products, Wilmington, DE). The results were further confirmed by analyzing the RNA on the Fragment Analyzer (Advanced Analytical Technologies, Ankeny, IA) where the samples showed RIN values ranging from 5 to 7 with only slight RNA degradation. Unfortunately, as this species is endangered and the government will only permit tissues to be harvested from animals that died accidentally, this is currently the best material that can be obtained.

Illumina TruSeq RNA library preparation and sequencing

Libraries were constructed following the manufacturer’s protocol with reagents supplied in Illumina’s TruSeq mRNA Stranded Library Preparation Kit. Briefly, the poly-A containing mRNA was purified from total RNA, mRNA was fragmented, double-stranded cDNA generated from fragmented RNA, and the index containing adapters were ligated to the ends.

The final construct of each purified library was evaluated using the Fragment Analyzer (Advanced Analytical Technologies, Ankeny, IA) automated electrophoresis system, quantified with the Qubit flourometer using the quant-iT HS dsDNA reagent kit (Invitrogen), and diluted according to Illumina’s standard sequencing protocol for sequencing on the NextSeq 500. Libraries were sequenced at the University of Missouri DNA Core Facility to obtain 75 base pair, single end reads.

Bioinformatics analyses

Since the Amami spiny rat genome has not been sequenced or annotated, we had to develop our own reference transcriptomic library. Previously described methods [97] were used to assemble a de novo reference transcriptome and analyze differential expression. Raw RNA-Seq reads were subjected to a multi-phase quality control regimen as previously described [98]. Specifically, 50-mer RNA-Seq reads from the Illumina HiSeq were first cleaned using the Fastx Toolkit (https://github.com/agordon/fastx_toolkit) to trim the 3′ ends of low quality (phred score < 20) bases and dropping reads when less than 32 bases remaining. Reads were then filtered to exclude those that did not have a minimum of 95% of their bases with a quality score of 20 or more. Adapter sequence was removed with CutAdapt, version 1.8.3 [99]. To create a final set of quality-controlled RNA-Seq reads, hereafter referred to as QC reads, foreign or undesirable reads were removed by sequence matching to the Phi-X genome (NC_001422.1) [100], the relevant ribosomal RNA genes as downloaded from the National Center for Biotechnology Information [101] or rodent repeat elements in RepBase, version 20.02 [102] using Bowtie, version 1.1.1 [103]. QC reads from all samples were used to generate a reference transcriptome assembly (RTA) by using Trinity, version 2.4.0, and default settings, except for random-access memory (RAM) and central processing unit (CPU) parameters [104]. Subsequently, QC reads from each sample were aligned to the RTA and expression estimates for each transcript > 200 nt in length were determined using RSEM via the align_and_estimate_abundance.pl script within the Trinity suite and these options: --prep_reference --est_method RSEM --aln_method bowtie2 --fragment_length 375 --fragment_std 40 --trinity_mode. The transcript expression estimates were normalized and prepared for downstream analysis using the abundance_estimates_to_matrix.pl script within the Trinity suite and these options: --est_method RSEM --cross_sample_norm TMM. Samples were subjected to a differential expression (DE) analysis using tximport [54] and DESeq2 [55]. Only transcripts having a sum of counts from all samples greater than 25 were analyzed. A transcript or gene was considered DE if the false discovery rate (FDR) associated with its expression ratio between conditions was less than 0.05. For enrichment, pathway and network analyses, we used various R packages and web sites, including EGSEA [105], DAVID [106, 107], TargetMine [108] and STRING [56]. Prior to using these programs, gene names were assigned to each DE transcript based on a BLASTX best hit strategy. Briefly, a gene name was assigned to a transcript if its best BLASTX alignment to either a human, rat or mouse protein carried an E-value of less than 1e− 06 and the aligned regions covered at least 40% of the reference protein. This protocol allowed us to assign gene names to approximately 65% of the differentially expressed transcripts. Assembly metrics were generated using TransRate [109] using the --assembly option and BUSCO, v3, using the --mode transcriptome option and the vertebrate database [49]. To analyze for transcript differences, we used the abundance_estimates_to_matrix.pl script included with the Trinity software to generate a matrix of counts for each transcript using the RSEM method. Subsequently, we used the run_DE_analysis.pl script included with the Trinity software to determine differentially expressed transcripts using the DESeq2 method. For each analysis, we used the list of differentially expressed (q < 0.05) gene symbols as input and with default options for the particular piece of software.

PCA plots were generated for genes and transcripts that showed sex dependent differences. The first and second principal components (which contain the two highest proportions of variation) were plotted on the X and Y axes, respectively. We used normalized gene expression estimates generated by the Trinity script abundance_estimates_to_matrix.pl using the TMM (Trimmed Mean of M values, [110]) normalization method.

qPCR analyses

Total RNA was reverse transcribed into cDNA using the QuantiTect Reverse Transcription Kit (Catalogue #205310, Qiagen). The quantitative polymerase chain reaction (qPCR) procedure was performed on the Applied Biosystems 7500 Real-Time PCR System (Carlsbad, CA) using the QuantiTect SYBR Green PCR Kit (Catalogue #204143; Qiagen). Primer sequences for the genes examined are listed in Additional file 6: Table S1, and primers were purchased from IDT (Coralville, IA). The qPCR conditions employed were 1) 15 min at 95 °C for polymerase activation 2) 40 cycles of: denaturation 40 s at 94 °C, annealing 40 s at 55 °C, and extension 72 °C for 1.50 min 3) Dissociation melt curve analysis from 60 °C to 90 °C. The internal control primer was glyceraldehyde-3-phosphate dehydrogenase (Gapdh) and test genes included: Svs5, Cyp1b1, and Serpina3n.

Statistics

The qPCR Data were analyzed by using GraphPad Prism 5 (GraphPad Software, La Jolla, CA). ANOVA was done to compare delta cycle threshold (dCt) values obtained in the brain samples of males to those obtained from female brain samples. However, in Fig. 6, data are expressed as relative fold change based on the 2-ΔΔCt method.

Abbreviations

AVPV:

Anteroventral periventricular nucleus

Bdkrb2:

B2 bradykinin receptor isoform XI

byr3:

Byr3 zinc finger protein

Cbln4:

Cerebellin 4 precursor gene

Cbx2:

Chromobox protein homolog 2

Col3a1:

Collagen type III alpha 1 chain

CPU:

Central processing unit

Csrp3:

Cysteine and glycine rich protein 3

Cul1:

Cullin 1

Cul3:

Cullin 3

Cyp1b1:

Cytochrome P450 1B1

Dcn:

Decorin

dCt:

Delta cycle threshold

DE:

Differential expression

Dmrt1:

Doublesex and mab-3 related transcription factor 1

E3:

Ubiquitin ligases

Eif2s3y:

Eukaryotic translation initiation factor 2 subunit 3, Y-linked

Er71:

Ets variant 2, Etv2

FCG:

Four core genotype

FDR:

False discovery rate

Fle1:

Flexuosa

Fmod:

Fibromodulin

GABA:

Gamma-aminobutyric acid

Gapdh:

Glyceraldehyde-3-phosphate dehydrogenase

GH:

Growth hormone

Hdac2:

Histone deacetylase 2

her:

Hermaphrodite

HP1:

Heterochromatin protein 1

iPSC:

Induced pluripotent stem cells

Kat6al:

Lysine acetyltransferase 6A-like

Kdm5d:

Lysine Demethylase 5D

KRAB:

Krueppel-associated box

KRAB-O:

Krueppel-associated box only

Lum:

Lumican precursor

Mb:

Myoglobin

Mbd4:

Methyl-CpG-binding domain protein 4

Mfn1:

Mitofusin 1

Mfn2:

Mitofusin 2

Mgrn1:

Mahogunin ring finger 1

Mki67:

Antigen KI-67 isoform X3

MPOA:

Medial preoptic area

Myl1:

Myosin light chain 1

Myl2:

Myosin regulatory light chain 2

Nt:

Nucleotides

Nt3:

Neurotrophin 3

Ogn:

Osteoglycin

Omd:

Osteomodulin

PCA:

Principal component analyses

Pja2:

Praja ring finger ubiquitin ligase 2

PR:

Progesterone receptor

QC:

Quality-controlled

Rag2:

V(D)J recombination-activating protein 2

RAM:

Random-access memory

Rnf220:

Ring finger protein 220

Rnf25:

Ring finger protein 25

RTA:

Reference transcriptome assembly

Sap30l:

Sin3-Associated Protein P30-Like

Serpina3n:

Serine (or cysteine) peptidase inhibitor, clade A, member 3 N

Setmar:

SET domain and mariner transposase fusion gene

Shprh:

SNF2 histone linker PHD RING helicase

Siah2:

Siah E3 Ubiquitin Protein Ligase 2

Smpx:

Small muscle protein X-linked

Sox9:

Sry-related box 9

Spata7:

Spermatogenesis Associated 7

SRY:

Sex-determining region Y

Svs5:

Seminal vesicle secretory protein 5

Tcf21:

Transcription factor 21

TDF:

Testis determining factor

Tex9:

Testis-expressed protein 9

Trdn:

Triadin

Zbtb44:

Zinc finger and BTB domain-containing protein 44

Zfp112:

Zinc finger protein 112

Zfp35:

Zinc finger protein 35

Zfp410:

Zinc finger protein 410

Zfp64:

Zinc finger protein 64

Zfp655:

Zinc finger protein 655

Zfp708:

Zinc finger protein 708

Zfx:

Zinc finger protein X-linked

Zfy:

Zinc finger protein Y-linked

Zpr1:

Zinc finger protein

References

  1. 1.

    Morris JA, Jordan CL, Breedlove SM. Sexual differentiation of the vertebrate nervous system. Nat Neurosci. 2004;7(10):1034–9.

  2. 2.

    Phoenix CH, Goy RW, Gerall AA, Young WC. Organizing action of prenatally administered testosterone propionate on the tissues mediating mating behavior in the female Guinea pig. Endocrinology. 1959;65(3):369–82.

  3. 3.

    Arnold AP, Breedlove SM. Organizational and activational effects of sex steroids on brain and behavior: a reanalysis. Horm Behav. 1985;19(4):469–98.

  4. 4.

    Watson J, Adkins-Regan E. Activation of sexual behavior by implantation of testosterone propionate and estradiol benzoate into the preoptic area of the male Japanese quail (Coturnix japonica). Horm Behav. 1989;23:251–68.

  5. 5.

    Watson J, Adkins-Regan E. Testosterone implanted in the preoptic area of male Japanese quail must be aromatized to activate copulation. Horm Behav. 1989;23:432–47.

  6. 6.

    Bowers JM, Waddell J, McCarthy MM. A developmental sex difference in hippocampal neurogenesis is mediated by endogenous oestradiol. Biol Sex Differ. 2010;1(1):8.

  7. 7.

    Konkle AT, McCarthy MM. Developmental time course of estradiol, testosterone, and dihydrotestosterone levels in discrete regions of male and female rat brain. Endocrinology. 2011;152(1):223–35.

  8. 8.

    Chung WC, Auger AP. Gender differences in neurodevelopment and epigenetics. Pflugers Arch. 2013;465(5):573–84.

  9. 9.

    Tsai HW, Grant PA, Rissman EF. Sex differences in histone modifications in the neonatal mouse brain. Epigenetics. 2009;4(1):47–53.

  10. 10.

    Jasarevic E, Geary DC, Rosenfeld CS. Sexually selected traits: a fundamental framework for studies on behavioral epigenetics. ILAR J. 2012;53(3–4):253–69.

  11. 11.

    Menger Y, Bettscheider M, Murgatroyd C, Spengler D. Sex differences in brain epigenetics. Epigenomics. 2010;2(6):807–21.

  12. 12.

    Matsuda KI, Mori H, Kawata M. Epigenetic mechanisms are involved in sexual differentiation of the brain. Rev Endocr Metab Disord. 2012;13(3):163–71.

  13. 13.

    Nugent BM, McCarthy MM. Epigenetic underpinnings of developmental sex differences in the brain. Neuroendocrinol. 2011;93(3):150–8.

  14. 14.

    Cortez D, Marin R, Toledo-Flores D, Froidevaux L, Liechti A, Waters PD, Grutzner F, Kaessmann H. Origins and functional evolution of Y chromosomes across mammals. Nature. 2014;508(7497):488–93.

  15. 15.

    Delbridge ML, Graves JA. Mammalian Y chromosome evolution and the male-specific functions of Y chromosome-borne genes. Rev Reprod. 1999;4(2):101–9.

  16. 16.

    Graves JA. Evolution of the testis-determining gene--the rise and fall of SRY. Novartis Found Symp. 2002;244:86–97 discussion 97–101, 203–106, 253–107.

  17. 17.

    Marshall Graves JA. The rise and fall of SRY. Trends Genet. 2002;18(5):259–64.

  18. 18.

    Graves JA. Weird mammals provide insights into the evolution of mammalian sex chromosomes and dosage compensation. J Genet. 2015;94(4):567–74.

  19. 19.

    Berta P, Hawkins JR, Sinclair AH, Taylor A, Griffiths BL, Goodfellow PN, Fellous M. Genetic evidence equating SRY and the testis-determining factor. Nature. 1990;348(6300):448–50.

  20. 20.

    Sinclair AH, Berta P, Palmer MS, Hawkins JR, Griffiths BL, Smith MJ, Foster JW, Frischauf AM, Lovell-Badge R, Goodfellow PN. A gene from the human sex- determining region encodes a protein with homology to a conserved DNA-binding motif. Nature. 1990;346(6281):240–4.

  21. 21.

    Foster JW, Brennan FE, Hampikian GK, Goodfellow PN, Sinclair AH, Lovell-Badge R, Selwood L, Renfree MB, Cooper DW, Graves JA. Evolution of sex determination and the Y chromosome: SRY-related sequences in marsupials. Nature. 1992;359(6395):531–3.

  22. 22.

    Hacker A, Capel B, Goodfellow P, Lovell-Badge R. Expression of Sry, the mouse sex determining gene. Development. 1995;121(6):1603–14.

  23. 23.

    Sekido R, Lovell-Badge R. Sex determination involves synergistic action of SRY and SF1 on a specific Sox9 enhancer. Nature. 2008;453(7197):930–4.

  24. 24.

    Bradford ST, Hiramatsu R, Maddugoda MP, Bernard P, Chaboissier MC, Sinclair A, Schedl A, Harley V, Kanai Y, Koopman P, et al. The cerebellin 4 precursor gene is a direct target of SRY and SOX9 in mice. Biol Reprod. 2009;80(6):1178–88.

  25. 25.

    Bhandari RK, Sadler-Riggleman I, Clement TM, Skinner MK. Basic helix-loop-helix transcription factor TCF21 is a downstream target of the male sex determining gene SRY. PLoS One. 2011;6(5):e19935.

  26. 26.

    Clement TM, Bhandari RK, Sadler-Riggleman I, Skinner MK. SRY directly regulates the neurotrophin 3 promoter during male sex determination and testis development in rats. Biol Reprod. 2011;85(2):277–84.

  27. 27.

    Wagner CK, Xu J, Pfau JL, Quadros PS, De Vries GJ, Arnold AP. Neonatal mice possessing an Sry transgene show a masculinized pattern of progesterone receptor expression in the brain independent of sex chromosome status. Endocrinology. 2004;145(3):1046–9.

  28. 28.

    Arnold AP, Chen X. What does the “four core genotypes” mouse model tell us about sex differences in the brain and other tissues? Front Neuroendocrinol. 2009;30(1):1–9.

  29. 29.

    Chen X, McClusky R, Chen J, Beaven SW, Tontonoz P, Arnold AP, Reue K. The number of x chromosomes causes sex differences in adiposity in mice. PLoS Genet. 2012;8(5):e1002709.

  30. 30.

    Kopsida E, Lynn PM, Humby T, Wilkinson LS, Davies W. Dissociable effects of Sry and sex chromosome complement on activity, feeding and anxiety-related behaviours in mice. PLoS One. 2013;8(8):e73699.

  31. 31.

    Kuljis DA, Loh DH, Truong D, Vosko AM, Ong ML, McClusky R, Arnold AP, Colwell CS. Gonadal- and sex-chromosome-dependent sex differences in the circadian system. Endocrinology. 2013;154(4):1501–12.

  32. 32.

    Gatewood JD, Wills A, Shetty S, Xu J, Arnold AP, Burgoyne PS, Rissman EF. Sex chromosome complement and gonadal sex influence aggressive and parental behaviors in mice. J Neurosci. 2006;26(8):2335–42.

  33. 33.

    Seney ML, Ekong KI, Ding Y, Tseng GC, Sibille E. Sex chromosome complement regulates expression of mood-related genes. Biol Sex Differ. 2013;4(1):20.

  34. 34.

    Quinnies KM, Bonthuis PJ, Harris EP, Shetty SR, Rissman EF. Neural growth hormone: regional regulation by estradiol and/or sex chromosome complement in male and female mice. Biol Sex Differ. 2015;6:8.

  35. 35.

    McPhie-Lalmansingh AA, Tejada LD, Weaver JL, Rissman EF. Sex chromosome complement affects social interactions in mice. Horm Behav. 2008;54(4):565–70.

  36. 36.

    Soullier S, Hanni C, Catzeflis F, Berta P, Laudet V. Male sex determination in the spiny rat Tokudaia osimensis (Rodentia: Muridae) is not Sry dependent. Mamm Genome. 1998;9(7):590–2.

  37. 37.

    Just W, Rau W, Vogel W, Akhverdian M, Fredga K, Graves JA, Lyapunova E. Absence of Sry in species of the vole Ellobius. Nat Genet. 1995;11(2):117–8.

  38. 38.

    Vogel W, Jainta S, Rau W, Geerkens C, Baumstark A, Correa-Cerro LS, Ebenhoch C, Just W. Sex determination in Ellobius lutescens: the story of an enigma. Cytogenet Cell Genet. 1998;80(1–4):214–21.

  39. 39.

    Matthey R. New type of chromosomal sex determination in the mammals Ellobius lutescens Th. And Microtus (Chilotus) oregoni Bachm. (Muridae, Microtinae). Experientia. 1958;14(7):240–1.

  40. 40.

    Fredga K. Aberrant chromosomal sex-determining mechanisms in mammals, with special reference to species with XY females. Philos Trans R Soc Lond Ser B Biol Sci. 1988;322(1208):83–95.

  41. 41.

    Sutou S, Mitsui Y, Tsuchiya K. Sex determination without the Y chromosome in two Japanese rodents Tokudaia osimensis osimensis and Tokudaia osimensis spp. Mamm Genome. 2001;12(1):17–21.

  42. 42.

    Kuroiwa A, Handa S, Nishiyama C, Chiba E, Yamada F, Abe S, Matsuda Y. Additional copies of CBX2 in the genomes of males of mammals lacking SRY, the Amami spiny rat (Tokudaia osimensis) and the Tokunoshima spiny rat (Tokudaia tokunoshimensis). Chromosom Res. 2011;19(5):635–44.

  43. 43.

    Otake T, Kuroiwa A. Molecular mechanism of male differentiation is conserved in the SRY-absent mammal, Tokudaia osimensis. Sci Rep. 2016;6:32874.

  44. 44.

    Gershoni M, Pietrokovski S. The landscape of sex-differential transcriptome and its consequent selection in human adults. BMC Biol. 2017;15(1):7.

  45. 45.

    Shi L, Zhang Z, Su B. Sex biased gene expression profiling of human brains at major developmental stages. Sci Rep. 2016;6:21181.

  46. 46.

    Trabzuni D, Ramasamy A, Imran S, Walker R, Smith C, Weale ME, Hardy J, Ryten M. Widespread sex differences in gene expression and splicing in the adult human brain. Nat Commun. 2013;4:2771.

  47. 47.

    Naurin S, Hansson B, Hasselquist D, Kim YH, Bensch S. The sex-biased brain: sexual dimorphism in gene expression in two species of songbirds. BMC Genomics. 2011;12:37.

  48. 48.

    Bundy JL, Vied C, Nowakowski RS. Sex differences in the molecular signature of the developing mouse hippocampus. BMC Genomics. 2017;18(1):237.

  49. 49.

    Simao FA, Waterhouse RM, Ioannidis P, Kriventseva EV, Zdobnov EM. BUSCO: assessing genome assembly and annotation completeness with single-copy orthologs. Bioinformatics. 2015;31(19):3210–2.

  50. 50.

    Khudyakov JI, Champagne CD, Meneghetti LM, Crocker DE. Blubber transcriptome response to acute stress axis activation involves transient changes in adipogenesis and lipolysis in a fasting-adapted marine mammal. Sci Rep. 2017;7:42110.

  51. 51.

    Mamrot J, Legaie R, Ellery SJ, Wilson T, Seemann T. De novo transcriptome assembly for the spiny mouse (Acomys cahirinus). bioRxiv. 2017;7(1):8996.

  52. 52.

    Morey JS, Burek Huntington KA, Campbell M, Clauss TM, Goertz CE, Hobbs RC, Lunardi D, Moors AJ, Neely MG, Schwacke LH, et al. De novo transcriptome assembly and RNA-Seq expression analysis in blood from beluga whales of Bristol Bay, AK. Mar Genomics. 2017;35:77–92.

  53. 53.

    MacManes MD. Severe acute dehydration in a desert rodent elicits a transcriptional response that effectively prevents kidney injury. Am J Physiol Renal Physiol. 2017;313(2):F262–72.

  54. 54.

    Soneson C, Love MI, Robinson MD. Differential analyses for RNA-seq: transcript-level estimates improve gene-level inferences. F1000Res. 2015;4:1521.

  55. 55.

    Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15(12):550.

  56. 56.

    Szklarczyk D, Morris JH, Cook H, Kuhn M, Wyder S, Simonovic M, Santos A, Doncheva NT, Roth A, Bork P, et al. The STRING database in 2017: quality-controlled protein-protein association networks, made broadly accessible. Nucleic Acids Res. 2017;45(D1):D362–8.

  57. 57.

    Zenteno JC, Lopez M, Vera C, Mendez JP, Kofman-Alfaro S. Two SRY-negative XX male brothers without genital ambiguity. Hum Genet. 1997;100(5–6):606–10.

  58. 58.

    Hyon C, Chantot-Bastaraud S, Harbuz R, Bhouri R, Perrot N, Peycelon M, Sibony M, Rojo S, Piguel X, Bilan F, et al. Refining the regulatory region upstream of SOX9 associated with 46,XX testicular disorders of sex development (DSD). Am J Med Genet A. 2015;167a(8):1851–8.

  59. 59.

    Vetro A, Ciccone R, Giorda R, Patricelli MG, Della Mina E, Forlino A, Zuffardi O. XX males SRY negative: a confirmed cause of infertility. J Med Genet. 2011;48(10):710–2.

  60. 60.

    Switonski M, Jackowiak H, Godynicki S, Klukowska J, Borsiak K, Urbaniak K. Familial occurrence of pig intersexes (38,XX; SRY-negative) on a commercial fattening farm. Anim Reprod Sci. 2002;69(1–2):117–24.

  61. 61.

    Rosenfeld CS. Brain sexual differentiation and requirement of SRY: why or why not? Front Neurosci. 2017;11:632.

  62. 62.

    Mayer A, Lahr G, Swaab DF, Pilgrim C, Reisert I. The Y-chromosomal genes SRY and ZFY are transcribed in adult human brain. Neurogenetics. 1998;1(4):281–8.

  63. 63.

    Clepet C, Schafer AJ, Sinclair AH, Palmer MS, Lovell-Badge R, Goodfellow PN. The human SRY transcript. Hum Mol Genet. 1993;2(12):2007–12.

  64. 64.

    Kawamata M, Inoue H, Nishimori K. Male-specific function of Dmrt7 by sexually dimorphic translation in mouse testis. Sex Dev. 2007;1(5):297–304.

  65. 65.

    Uto K, Nakajo N, Sagata N. Two structural variants of Nek2 kinase, termed Nek2A and Nek2B, are differentially expressed in Xenopus tissues and development. Dev Biol. 1999;208(2):456–64.

  66. 66.

    Cho S, Huang ZY, Zhang J. Sex-specific splicing of the honeybee doublesex gene reveals 300 million years of evolution at the bottom of the insect sex-determination pathway. Genetics. 2007;177(3):1733–41.

  67. 67.

    Duan J, Xu H, Guo H, O'Brochta DA, Wang F, Ma S, Zhang L, Zha X, Zhao P, Xia Q. New insights into the genomic organization and splicing of the doublesex gene, a terminal regulator of sexual differentiation in the silkworm Bombyx mori. PLoS One. 2013;8(11):e79703.

  68. 68.

    Oh HJ, Li Y, Lau YF. Sry associates with the heterochromatin protein 1 complex by interacting with a KRAB domain protein. Biol Reprod. 2005;72(2):407–15.

  69. 69.

    Oh HJ, Lau YF. KRAB: a partner for SRY action on chromatin. Mol Cell Endocrinol. 2006;247(1–2):47–52.

  70. 70.

    Su H, Lau YF. Demonstration of a stage-specific expression of the ZFY protein in fetal mouse testis using anti-peptide antibodies. Mol Reprod Dev. 1992;33(3):252–8.

  71. 71.

    Lopez M, Torres L, Mendez JP, Cervantes A, Alfaro G, Perez-Palacios G, Erickson RP, Kofman-Alfaro S. SRY alone can induce normal male sexual differentiation. Am J Med Genet. 1995;55(3):356–8.

  72. 72.

    Kuroiwa A, Ishiguchi Y, Yamada F, Shintaro A, Matsuda Y. The process of a Y-loss event in an XO/XO mammal, the Ryukyu spiny rat. Chromosoma. 2010;119(5):519–26.

  73. 73.

    Luoh SW, Bain PA, Polakiewicz RD, Goodheart ML, Gardner H, Jaenisch R, Page DC. Zfx mutation results in small animal size and reduced germ cell number in male and female mice. Development. 1997;124(11):2275–84.

  74. 74.

    Li H, Baker BS. Her, a gene required for sexual differentiation in Drosophila, encodes a zinc finger protein with characteristics of ZFY-like proteins and is expressed independently of the sex determination hierarchy. Development. 1998;125(2):225–35.

  75. 75.

    Coppin E. The fle1 gene encoding a C2H2 zinc finger protein co-ordinates male and female sexual differentiation in Podospora anserina. Mol Microbiol. 2002;43(5):1255–68.

  76. 76.

    Xu HP, Rajavashisth T, Grewal N, Jung V, Riggs M, Rodgers L, Wigler M. A gene encoding a protein with seven zinc finger domains acts on the sexual differentiation pathways of Schizosaccharomyces pombe. Mol Biol Cell. 1992;3(7):721–34.

  77. 77.

    Okada G, Maruo K, Funada S, Nakamura M. Differential display analysis of gene expression in female-to-male sex-reversing gonads of the frog Rana rugosa. Gen Comp Endocrinol. 2008;155(3):623–34.

  78. 78.

    Cunliffe V, Williams S, Trowsdale J. Genomic analysis of a mouse zinc finger gene, Zfp-35, that is up-regulated during spermatogenesis. Genomics. 1990;8(2):331–9.

  79. 79.

    Cunliffe V, Koopman P, McLaren A, Trowsdale J. A mouse zinc finger gene which is transiently expressed during spermatogenesis. EMBO J. 1990;9(1):197–205.

  80. 80.

    Sato Y, Shinka T, Chen G, Yan HT, Sakamoto K, Ewis AA, Aburatani H, Nakahori Y. Proteomics and transcriptome approaches to investigate the mechanism of human sex determination. Cell Biol Int. 2009;33(8):839–47.

  81. 81.

    Gelfman S, Ast G. When epigenetics meets alternative splicing: the roles of DNA methylation and GC architecture. Epigenomics. 2013;5(4):351–3.

  82. 82.

    Gelfman S, Cohen N, Yearim A, Ast G. DNA-methylation effect on cotranscriptional splicing is dependent on GC architecture of the exon-intron structure. Genome Res. 2013;23(5):789–99.

  83. 83.

    Luco RF, Allo M, Schor IE, Kornblihtt AR, Misteli T. Epigenetics in alternative pre-mRNA splicing. Cell. 2011;144(1):16–26.

  84. 84.

    Maunakea AK, Chepelev I, Cui K, Zhao K. Intragenic DNA methylation modulates alternative splicing by recruiting MeCP2 to promote exon recognition. Cell Res. 2013;23(11):1256–69.

  85. 85.

    Schor IE, Fiszbein A, Petrillo E, Kornblihtt AR. Intragenic epigenetic changes modulate NCAM alternative splicing in neuronal differentiation. EMBO J. 2013;32(16):2264–74.

  86. 86.

    Yearim A, Gelfman S, Shayevitch R, Melcer S, Glaich O, Mallm JP, Nissim-Rafinia M, Cohen AH, Rippe K, Meshorer E, et al. HP1 is involved in regulating the global impact of DNA methylation on alternative splicing. Cell Rep. 2015;10(7):1122–34.

  87. 87.

    Ge C, Ye J, Zhang H, Zhang Y, Sun W, Sang Y, Capel B, Qian G. Dmrt1 induces the male pathway in a turtle species with temperature-dependent sex determination. Development. 2017;144(12):2222–33.

  88. 88.

    Matsumoto Y, Hannigan B, Crews D. Embryonic PCB exposure alters phenotypic, genetic, and epigenetic profiles in turtle sex determination, a biomarker of environmental contamination. Endocrinology. 2014;155(11):4168–77.

  89. 89.

    Matsumoto Y, Hannigan B, Crews D. Temperature shift alters DNA methylation and histone modification patterns in gonadal aromatase (cyp19a1) gene in species with temperature-dependent sex determination. PLoS One. 2016;11(11):e0167362.

  90. 90.

    Radhakrishnan S, Literman R, Mizoguchi B, Valenzuela N. MeDIP-seq and nCpG analyses illuminate sexually dimorphic methylation of gonadal development genes with high historic methylation in turtle hatchlings with temperature-dependent sex determination. Epigenetics Chromatin. 2017;10:28.

  91. 91.

    Venegas D, Marmolejo-Valencia A, Valdes-Quezada C, Govenzensky T, Recillas-Targa F, Merchant-Larios H. Dimorphic DNA methylation during temperature-dependent sex determination in the sea turtle Lepidochelys olivacea. Gen Comp Endocrinol. 2016;236:35–41.

  92. 92.

    Honda A. Applying iPSCs for preserving endangered species and elucidating the evolution of mammalian sex determination. Bioessays. 2018;40(6):e1700152.

  93. 93.

    Honda A, Choijookhuu N, Izu H, Kawano Y, Inokuchi M, Honsho K, Lee AR, Nabekura H, Ohta H, Tsukiyama T, et al. Flexible adaptation of male germ cells from female iPSCs of endangered Tokudaia osimensis. Science Adv. 2017;3(5):e1602179.

  94. 94.

    Murata C, Kuroki Y, Imoto I, Kuroiwa A. Ancestral Y-linked genes were maintained by translocation to the X and Y chromosomes fused to an autosomal pair in the Okinawa spiny rat Tokudaia muenninki. Chromosom Res. 2016;24(3):407–9.

  95. 95.

    Murata C, Yamada F, Kawauchi N, Matsuda Y, Kuroiwa A. Multiple copies of SRY on the large Y chromosome of the Okinawa spiny rat, Tokudaia muenninki. Chromosome Res. 2010;18(6):623–34.

  96. 96.

    Murata C, Yamada F, Kawauchi N, Matsuda Y, Kuroiwa A. The Y chromosome of the Okinawa spiny rat, Tokudaia muenninki, was rescued through fusion with an autosome. Chromosom Res. 2012;20(1):111–25.

  97. 97.

    Johnson SA, Spollen WG, Manshack LK, Bivens NJ, Givan SA, Rosenfeld CS. Hypothalamic transcriptomic alterations in male and female California mice (Peromyscus californicus) developmentally exposed to bisphenol A or ethinyl estradiol. Physiol Rep. 2017;5(3). https://doi.org/10.14814/phy2.13133.

  98. 98.

    Givan SA, Bottoms CA, Spollen WG. Computational analysis of RNA-seq. Methods Mol Biol. 2012;883:201–19.

  99. 99.

    Martin M. Cutadapt removes adaptor sequences from high-throughput sequencing reads. EMBnetjournal. 2011;17:10.

  100. 100.

    Eppsl-I-G-N: Available from: http://www.ncbi.nlm.nih.gov/genome/?term=NC_001422.1.

  101. 101.

    NCfBI: Available from: https://www.ncbi.nlm.nih.gov/.

  102. 102.

    Jurka J, Kapitonov VV, Pavlicek A, Klonowski P, Kohany O, Walichiewicz J. Repbase update, a database of eukaryotic repetitive elements. Cytogenet Genome Res. 2005;110(1–4):462–7.

  103. 103.

    Langmead B, Trapnell C, Pop M, Salzberg SL. Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biol. 2009;10(3):R25.

  104. 104.

    Haas BJ, Papanicolaou A, Yassour M, Grabherr M, Blood PD, Bowden J, Couger MB, Eccles D, Li B, Lieber M, et al. De novo transcript sequence reconstruction from RNA-seq using the trinity platform for reference generation and analysis. Nat Protoc. 2013;8(8):1494–512.

  105. 105.

    Alhamdoosh M, Ng M, Wilson NJ, Sheridan JM, Huynh H, Wilson MJ, et al. Combining multiple tools outperforms individual methods in gene set enrichment analyses. bioRxiv. 2016;11:042580.

  106. 106.

    Huang DW, Sherman BT, Lempicki RA. Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat Protoc. 2008;4(1):44–57.

  107. 107.

    Huang da W, Sherman BT, Lempicki RA. Bioinformatics enrichment tools: paths toward the comprehensive functional analysis of large gene lists. Nucleic Acids Res. 2009;37(1):1–13.

  108. 108.

    Chen YA, Tripathi LP, Mizuguchi K. An integrative data analysis platform for gene set analysis and knowledge discovery in a data warehouse framework. Database. 2016;2016:baw009.

  109. 109.

    Smith-Unna R, Boursnell C, Patro R, Hibberd JM, Kelly S. TransRate: reference-free quality assessment of de novo transcriptome assemblies. Genome Res. 2016;26(8):1134–44.

  110. 110.

    Robinson MD, Oshlack A. A scaling normalization method for differential expression analysis of RNA-seq data. Genome Biol. 2010;11(3):R25.

Download references

Acknowledgements

Gratitude is extended to all those involved in collecting brain samples from Amami spiny rats and Sarah A. Johnson for assisting in the RNA isolation.

Funding

Support was provided by the University of Missouri Informatics Research Core Facility and Bond Life Sciences Center.

Availability of data and materials

The datasets supporting the conclusions of this article are included within the article (and its additional files). The reference transcriptomic library generated for Amami spiny rat as part of these studies is available at: https://www.ncbi.nlm.nih.gov/bioproject/474959.

Author information

TJ and AK obtained permissions to collect the samples and shipped them to the University of Missouri. MTO and NJB isolated the RNA and performed the RNAseq analyses. SAG analyzed the data. The manuscripts was written and revised by MTO, NJB, TK, AK, SAG, and CSR. All authors have read and approved the final manuscript.

Correspondence to Cheryl S. Rosenfeld.

Ethics declarations

Ethics approval and consent to participate

T. osimensis is endangered (The IUCN Red List of Threatened Species; https://www.iucnredlist.org/). This species has been protected by the Japanese government as a natural monument since 1972 and Nationally Endangered Species of Wild Fauna and Flora since 2016. With permission from the Agency for Cultural Affairs and the Ministry of the Environment in Japan, tissues were harvested from animals died accidentally during a mongoose eradication project on the island in which this species inhabits. In such cases of accidental death, the Japanese government has approved harvesting of brain and other tissues. A copy of this governmental approval, which has also been translated into English, has been provided to the journal.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Additional files

Additional file 1:

Gene expression difference in the brain of male vs. female Amami spiny rats. (XLSX 29 kb)

Additional file 2:

Transcript isoforms differentially regulated in the brain of male vs. female Amami spiny rats. (XLSX 529 kb)

Additional file 3:

Pathways upregulated in male Amami spiny rats based on differential gene expression. (XLSX 10 kb)

Additional file 4:

Pathways upregulated in male Amami spiny rats based on differential transcript isoform expression. (XLSX 30 kb)

Additional file 5:

Pathways upregulated in female Amami spiny rats based on differential transcript isoform expression. (XLSX 62 kb)

Additional file 6:

Table S1. Primer sequences used for gene expression analyses. (DOCX 14 kb)

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Keywords

  • Sexual differentiation
  • Organizational-Activational programming
  • Rodent
  • Endangered animals
  • SRY
  • Steroid hormones
  • Testosterone
  • Estrogen
  • RNAseq
  • Bioinformatics