Stress-mediated convergence of splicing landscapes in male and female rock doves
BMC Genomics volume 21, Article number: 251 (2020)
The process of alternative splicing provides a unique mechanism by which eukaryotes are able to produce numerous protein products from the same gene. Heightened variability in the proteome has been thought to potentiate increased behavioral complexity and response flexibility to environmental stimuli, thus contributing to more refined traits on which natural and sexual selection can act. While it has been long known that various forms of environmental stress can negatively affect sexual behavior and reproduction, we know little of how stress can affect the alternative splicing associated with these events, and less still about how splicing may differ between sexes. Using the model of the rock dove (Columba livia), our team previously uncovered sexual dimorphism in the basal and stress-responsive gene transcription of a biological system necessary for facilitating sexual behavior and reproduction, the hypothalamic-pituitary-gonadal (HPG) axis. In this study, we delve further into understanding the mechanistic underpinnings of how changes in the environment can affect reproduction by testing the alternative splicing response of the HPG axis to an external stressor in both sexes.
This study reveals dramatic baseline differences in HPG alternative splicing between males and females. However, after subjecting subjects to a restraint stress paradigm, we found a significant reduction in these differences between the sexes. In both stress and control treatments, we identified a higher incidence of splicing activity in the pituitary in both sexes as compared to other tissues. Of these splicing events, the core exon event is the most abundant form of splicing and more frequently occurs in the coding regions of the gene. Overall, we observed less splicing activity in the 3’UTR (untranslated region) end of transcripts than the 5’UTR or coding regions.
Our results provide vital new insight into sex-specific aspects of the stress response on the HPG axis at an unprecedented proximate level. Males and females uniquely respond to stress, yet exhibit splicing patterns suggesting a convergent, optimal splicing landscape for stress response. This information has the potential to inform evolutionary theory as well as the development of highly-specific drug targets for stress-induced reproductive dysfunction.
Organismal behavior and its mechanistic underpinnings have been consistent quandaries for many biologists [1,2,3,4]. Even among the few cases in which an adaptive behavior is clearly understood, we have little information on how proximate mechanisms have led to ultimate behavioral adaptations. One of the mechanisms by which animals modulate their response to stimuli is by adjusting the levels of endogenous protein products in physiological systems [5, 6]. This variable production may result in differential binding of hormones and signaling peptides, enabling an organism to receive information more accurately about external stimuli and effectively react [7, 8]. In addition to adjusting the quantity of a given protein, organisms may further alter these proteins by producing slightly modified versions (e.g., isoforms) of each transcript [9, 10]. This variable modification of gene products is called alternative splicing. Alternative splicing is a mechanism by which organisms can respond to their surroundings with extreme precision. Responding to stress requires this level of precision. As such, one can anticipate finding alternative splicing contributing to the organismal stress response.
Alternative splicing is a process common to eukaryotes  that involves cleavage of transcribed RNA at specific splice sites and varying inclusion or exclusion of genomic elements (introns and exons). In the human genome, approximately 80% of exons are > 200 bp in length [12, 13]; however, exon sizes identified in other species vary from a single base to > 17,000 bp in length. Each human gene contains, on average, eight exons . This variable inclusion of genetic sequences results in a dramatic increase in the number of potential transcript and protein products that a single gene may produce. Alternative splicing presents an additional mechanism by which mRNA levels and gene expression can be regulated, while also greatly increasing proteome diversity. Splicing activity is thought to be responsible for the majority of proteomic diversity in eukaryotes  and, potentially, may be an underlying mechanism of functional genomic evolution .
Numerous types of splicing events exist that occur at different frequencies in a given genome and alter proteins in subtle to dramatically different ways . Cassette exon splicing (also referred to as exon skipping) is the most common type of splicing event in vertebrates and invertebrates, while intron retention is more common in plants. Additionally, alternative selection of 5′ and 3′ splice sites, coupled with variable adenylation of the transcript, results in further modification of protein products [18, 19]. The splicing process consists of two major steps: assembly of the spliceosome and the actual splicing of pre-mRNA . In brief, the spliceosome is comprised of several small nuclear ribonucleoproteins that positionally establish the 5′ splice site, the branch point sequence, and the 3′ site. An assembly of spliceosome complexes and eight evolutionarily-conserved RNA-dependent (Ribonucleic Acid) ATPases/helicases (Adenosine Triphosphate) is then followed by the execution of numerous splicing steps, ultimately resulting in exon excision, exon ligation, or intron retention . The inclusion of an exon in the final mRNA product is entirely driven by cis- and trans-acting elements/factors. The interaction of these elements within the splicing process promotes or inhibits spliceosome activity on various splice regions, resulting in alternative splicing [21, 22].
Alternative splicing mechanisms enable organisms to sense and react to minute changes in the local environment, allowing both plants and animals to tailor their responses to their surroundings with extreme precision [23, 24]. Previous research has revealed unique roles for alternative splicing in the immune response of chickens with avian pathogenic E.coli , mediation of abiotic stress response pathways of plants , and enhanced fear memory of mice . Alternative splicing has also been implicated in various aspects of cancer, including oncogenesis  and cancer drug resistance [29, 30]. Some studies have identified a sex-bias in alternative splicing in Drosophila [31,32,33], while others have identified unique sex-specific splicing differences in human brains . The diverse roles of alternative splicing in biological processes and behavioral responses inherently speak to the depth and breadth that alternative splicing drives organismal physiology and behavior, at both local and global levels. By identifying the splicing landscape that modulates gene expression and mRNA transcript composition in both males and females, we increase the resolution at which we can comprehend the proximate mechanisms underlying animal physiology and behavior.
In vertebrates, a symphony of physiological events is required to regulate sexual behavior and reproduction, and these mechanisms are driven by an interconnected biological system made up of the hypothalamus in the brain, the pituitary gland, and the gonads (testes/ovaries) [7, 35,36,37]. This hypothalamic-pituitary-gonadal (HPG) axis can be disrupted in multiple, complex ways [7, 38,39,40]. However, we know little about how stress affects the HPG axis at the level of alternative splicing, and we know even less regarding its effects at this level in males versus females. Understanding how the alternative splicing landscape of the reproductive axis changes in the face of stress will not only offer more insight into how stress can affect reproduction, but deepen our proximate knowledge of biological processes and sexaully-biased behavioral responses in general.
Using the classic reproductive model [41,42,43,44] and rising genomics model [7, 35, 45,46,47] of the rock dove, Columba livia, we have identified sexual dimorphism in both basal  and restraint stress-responsive  HPG gene expression at the level of RNA transcription. In this study, we traverse beyond the level of transcription to test for sex-biased alternative splicing patterns in the HPG axis of the rock dove in response to a restraint stress stimulus. Using a relatively highly-replicated (n = 12/sex) study design, we identify significantly similar and different splicing events between the sexes and in response to restraint stress treatment. To our knowledge, this is the first report of sex-specific splicing events in the HPG axis in response to a stressor.
Sequencing results, read data, and code availability
Samples were sequenced to a read depth between 2.3 million and 24.5 million read pairs, for a total of 1,095,954,918 paired-end reads (more fully described in . Read data corresponding to the control birds are available using the European Nucleotide Archive project ID PRJEB16136; read data corresponding to the stressed birds are available at PRJEB21082. Read abundance and data on reads mapped can be found in Table S1. Additionally, genome annotation statistics can be found in Table S2. All code for analyses in this manuscript can be found at https://github.com/AndrewLangvt/Scripts/tree/master/splicing_analysis/.
Male vs. female splicing comparison
Our first aim was to understand sex-typical splicing in the hypothalamus (hyp) and pituitary (pit) by assessing each tissue for alternative splicing events between males and females. We counted the number, and type, of alternatively spliced loci between males and females in each treatment state (control: male vs female; stress: male vs female). This approach allowed us to determine how the splicing landscape changed between sexes in response to restraint stress, and in which state the sexes shared a more similar splicing profile. As previously stated, we did not include gonads in this comparison due to inherent splicing differences between tissue types. Chi-squared tests (hereafter, ChiSq) were used to determine statistical significance (p < 0.05) throughout our analyses; all p-values, degrees of freedom, and sample sizes are included in Table 1. Chi-squared tests were used to test null hypotheses that AS event abundance did not differ between treatments, sexes, tissue, type, or region (i.e that splicing events would be evenly distributed across whichever parameters we were considering).
Male vs. female splicing comparison: events by type
In total, we identified 158 splicing events in the hypothalamus and 225 events in the pituitary. When compared to the hypothalamus, the 42% increase of splicing event abundance seen in the pituitary is significant (ChiSq p = 6.18e-4). In both tissues, more events were identified in the control state compared to the stress condition (hyp: 99 control/59 stress, pit: 123 control/102 stress), but only the relationship in the hypothalamus was statistically significant (Fig. 1, ChiSq p: hyp = 1.46e-3; pit = 0.162).
These total counts were further broken down by event type (Fig. 1). The core exon event was the most abundant event identified across sex in both the hypothalamus and pituitary, regardless of treatment (ChiSq p: hyp = 2.22e-13, pit = 5.70e-43). The core exon event called by the software package. Whippet is a splicing event involving a full exonic segment: previously referred to as “cassette exon” or “exon skipping” in other publications. Of these core exon splice events, there were almost twice as many in the pituitary compared to hypothalamus (pit: 74 control/67 stress, hyp: 48 control/30 stress). Within these event types, we tested for statistical significance between splicing differences in males and females of each treatment group. Retained intron events in the hypothalamus were the only event to differ significantly between treatments (ChiSq p = 0.012), with nearly 3 times (280% increase) more splicing events in the control state than the stressed. Both core exon events in the hypothalamus and retained intron events in the pituitary reflected a similar increased abundance in splicing events of the control state, though these relationships were not significant (ChiSq p: CE-Hyp = 0.053, RI-Pit = 0.052). The distribution of Percent Spliced In (PSI) values between males and females did not vary between treatments, indicating that the level of event inclusion/exclusion difference between the sexes was generally unaffected by treatment.
Male vs. female splicing comparison: genes of interest
Using our comparison of male to female splicing patterns, we were able to identify sex-specific alternatively spliced genes in the stress response. We provide a full list of spliced genes within each comparison (Table S3) and also a complete list of all events including dPSI (delta PSI), probability, and genomic location (Table S4). Some of these spliced genes are involved in functional gene expression within the HPG axis. POU class 2 homeobox 1 (POU2F1), a transcription factor that regulates transcription of gonadotropin-releasing hormone (GnRH) [48, 49], is alternatively spliced in the male pituitary stress response. GnRH is a primary regulator of the HPG axis [50,51,52]. Splicing of POU2F1 likely affects the HPG axis, indirectly, by modulating transcription of GnRH [53,54,55]. Through alternative splicing of the POU2F1 gene in the pituitary, males may be altering signaling pathways within the HPG axis to optimize stress response.
Few genes were consistently alternatively spliced between males and females in both treatments. Those that did exhibit consistent alternative splicing between the sexes were often related to immune function. Rap1 GTP-ase activating protein 1 (RAP1GAP) is consistently alternatively spliced between male and female hypothalami, in both the control and stress treatments. Previous findings have shown RAP1GAP to be a putative oncogene . This gene mediates the strength of cell adhesions through regulation of Rap1, thus modulating T-cell response . Alternatively spliced between sexes in the pituitary, P-selectin (SELP) is known to preserve immune function in mice . The corresponding ligand, P-selectin glycoprotein ligand-1 (PSGL-1), negatively regulates T-cell response through binding of SELP . Genes consistently alternatively spliced between males and females may reflect splicing-level sexual dimorphism, indicating that males and females inherently differ in their splicing landscapes. Further, these differences appear to speak to a sex-specific stress response through modification of genes related to immune processes. Through future study of the splicing landscape of these genes of interest in additional tissues and states, we will likely reveal additional inherent splicing differences between the sexes.
Male vs. female splicing comparison: gene ontology
By observing abundances of parent ontology terms significantly deviating from genomic expectation, we were able to gain a broader understanding of gene-types targeted by alternative splicing. In the list of significant Molecular Function terms (Fig. 2), splicing of organic and heterocyclic compound genes is underrepresented in the pituitary, while small molecular and drug binding is overrepresented in the hypothalamus regardless of treatment. There was an interaction seen in the hypothalamus; splicing of organic and heterocyclic compound genes was overrepresented in the stress treatment, though underrepresented in the control group. Biological Process terms suggest that males and females differ very little in their splicing profile of metabolic genes in either the hypothalamus or pituitary, given there were fewer spliced genes with metabolism GO terms in these tissues (Figure S1). Finally, splicing events in stressed males and females are more abundant in genes related to cell/neuronal structure of the pituitary than the hypothalamus (Figure S2).
Male vs. female splicing comparison: exon size, location, and motifs
To further characterize the splicing profile of the sexes, we visualized distributions of exon sizes, where these exons were located, and protein motifs contained therein. In terms of size distribution, our analyses revealed no significant difference between the control (male vs female) and stress (male vs female) comparisons in the hypothalamus or pituitary, suggesting that the size of alternatively spliced exons between males and females does not differ between control or stressed states (Fig. 3). In all cases, spliced exons were smaller than predicted by the genomic distribution (Wilcoxon test p: hyp-control = 4.52e-12; hyp-stress = 2.16e-3; pit-control = 2.75e-8; pit-stress = 3.53e-10) (Fig. 3). Splicing events occurred in a variety of protein motifs, but no particular motif was significantly more spliced than the others (Figure S3).
In both tissues, core exon splice sites occurred in protein coding sequences much more frequently than either of the untranslated regions (ChiSq p: hyp = 5.95e-43, pit = 6.52e-141). This, perhaps, is not surprising given that the CDS regions are more abundant in the genome and alteration to these regions will ultimately result in changes to the protein sequence. The abundance of alternatively spliced exons present in the 5′ & 3′ untranslated regions of the hypothalamus and pituitary was significantly different from expected values (Fig. 4). In the hypothalamus, more spliced exons occurred in the 5’UTR than genomic proportions would predict, with more dramatic shifts in the restraint stress treatment than the control (ChiSq, p = 0.015). In the pituitary, the control group exhibited more spliced exons of the 5’UTR and both treatment groups presented more than a 6% decrease of splicing events in 3’UTR regions than predicted from genomic values (ChiSq p: 5’UTR = 0.041, 3’UTR = 2.60e-7).
Control vs. stress splicing comparison
The second aim of this study was to observe splicing differences within each sex in control and stress states, to identify how each sex individually responded to restraint stress. We compared, within each sex, control and stress treatments (female: control to stress; male: control to stress). We counted the number of alternative splicing events for these within-sex, across-treatment comparisons. Here, we include comparisons between gonads as the alternative splicing events identified are within-sex and thus, we can observe how male and female gonads individually respond to restraint stress. ChiSq tests were used to determine statistical significance throughout our analyses; all p-values, degrees of freedom, and sample sizes are included in Table 1.
Control vs. stress splicing comparison: events by type
Similar to our comparison of male-female alternative splicing, control-stress splicing reveals more splicing in the pituitary than other tissues (ChiSq p = 1.04e-6). The male hypothalamus is more active in stress-response splicing compared to the female hypothalamus (ChiSq p = 9.39e-4), and the female gonad exhibits more activity than the male gonad (ChiSq p = 1.92e-4). The male hypothalamus displays 59% more splicing events than the female hypothalamus, and we identified 78% more splicing events in the ovaries than the testes (Fig. 1).
Of all event types, the core exon event was most abundant in all tissues (ChiSq p: hyp = 3.64e-7, pit = 1.34e-43, gon = 8.74e-15) (Fig. 1). Of these core exon events, the gonads were the only tissue to present statistical significance across sexes; more core exon splice events occurred in the ovaries than in the testes (ChiSq p = 3.01e-3). This parallels our previous findings of elevated female gonadal gene expression in response to restraint stress . The only other event that differed significantly between the sexes was hypothalamic retained intron events, with more found in males than females (ChiSq p = 0.041).
Control vs. stress splicing comparison: genes of interest
Through assessment of alternative splicing events between control and stress states, we were able to identify genes spliced in response to stress within each sex. Estradiol 17-beta-dehydrogenase 11 (HSD17B11) was alternatively spliced between treatments in the male hypothalamus. HSD17B11 plays a role in hormone metabolism, through which it may mediate endogenous estrogen levels . Through feedback on the brain, estrogens control the pulsatile release of GnRH and can influence stress signaling  and enhance hypothalamic-pituitary-adrenal (HPA) function . HSD17B11 is also related to CEBPB (CCAAT/enhancer-binding protein beta), a transcription factor regulating the expression of genes involved in immune and inflammatory responses [63, 64]. The splicing of this gene, and others in our list, suggest a sex-specific response to stress mediated via alternative splicing.
As in our male-female splicing analysis, we identified few genes consistently alternatively spliced; however, the genes that were consistently spliced in both sexes were quite interesting. Rho family GTPases function as molecular switches to regulate a variety of cellular processes including cytoskeletal growth, gene transcription, and cytokinesis [65, 66]. Both males and females displayed alternative splicing of the pituitary gene Rho GTPase Activating Protein 32 (ARHGAP32). Through response to GnRH pulsatile expression, Rho GTPases regulates cell motility and cytoskeletal rearrangements [67, 68]. The shift from a control to stressed state likely places different constraints on cellular motility, requiring expedited movement of signaling molecules and energy transport throughout an organism; consistent splicing of ARHGAP32 may be one means by which cellular movement is adjusted during the stress response.
The hypothalamus of both males and females exhibit alternative splicing of the Retinitis Pigmentosa GTPase Regulator Interacting Protein 1 Like (RPGRIP1L) gene. RPGRIP1L is an adapter protein, localized at the ciliary transition zone, that interacts with G-protein coupled thromboxane A2 receptor (TBXA2R) to negatively regulate signaling at the ciliary base [69, 70]. This ciliary protein governs both autophagy and proteasome activity [69, 71]. As both sexes respond to stress by alternatively splicing the RPGRIP1L gene, the alteration of this locus may enable organisms to respond more effectively to a stressful environment through signal transduction and processing/disposal of damaged cells. Additional observation of the functioning of these gene products in control and stress states will provide further understanding as to the role these loci play in stress response.
Control vs. stress splicing comparison: gene ontology
In addition to characterizing splicing on a gene-by-gene basis, we assessed the gene lists, as a whole, to determine if a particular type of gene or genes related to a specific biological process were more readily spliced in the stress response. Our ontology analysis of the stress vs. control splicing events revealed elevated levels of splicing in genes with the following GO terms: small molecule binding, carb derivative binding, and drug binding (Fig. 5). These terms were more abundant than genomic predictions in the male hypothalamus, female pituitary, and male gonad, while these same terms were less abundant than expected in the female hypothalamus and female gonad. Of the male tissues, the hypothalamus and testes appear more active in this area. There also were more spliced genes related to drug binding than expected. Similar analyses were observed for Biological Function (Figure S4) and Cellular Component (Figure S5), highlighting sexual dimorphism in AS of hypothalamic and pituitary genes.
Control vs. stress splicing comparison: exon size, location, and motifs
Analysis of distributions of alternatively spliced exon sizes reveal that all tissues, regardless of sex, exhibit distributions of control-stress splicing events smaller than the genomic distribution (Wilcoxon, all p < 2.80e-3) (Fig. 3). The pituitary was the only tissue in which male and female distributions differed, with the male pituitary exons being longer than female (Wilcoxon, p = 0.022). Similar to the male-female splicing analysis of splicing in protein motifs, we found no significant difference in the distribution of events occurring in any one motif (Figure S6).
In analyzing locations of spliced exons, we again identified more events in coding regions than either 5′ or 3′ UTR, for all tissues (ChiSq p: hyp = 4.41e-40, pit = 9.60e-233, gon = 5.13e-67) (Fig. 6). Similar to our male-female analysis, we also identified a significantly lower abundance of splicing events in the 3’UTR regions of all tissues in the control-stress analysis (ChiSq p: hyp = 0.042, pit = 2.87e-12, gon = 8.77e-4). There were fewer events in the 5’UTR of the pituitary (ChiSq p = 2.04e-4), and more events occurred in the pituitary CDS than expected (ChiSq p = 3.50e-4).
Treatment-specific and sex-specific isoforms
Because we know that treatment and sex affect gene expression , we assessed each tissue for treatment-specific, and sex-specific isoforms of a gene. We did not identify any treatment-, nor sex-specific isoforms. This result suggests that even though gene expression may shift, and splicing events can be specific to an external stimulus, specific isoforms are not exclusive to one sex or treatment. A potential shortcoming of this particular analysis is the heavy reliance upon annotation. As the C. livia genome continues to improve with the progression of technology and sequence analysis tools we encourage others to reassess this work.
Bimodal distribution of PSI
We assessed the distributions of PSI values for each splicing event type across all analysis comparisons to determine if there existed a differential in the frequency that a splicing event was present within a given state. The measure of “percent spliced in” reflects the frequency that one splicing isoform was seen over another in an individual. Ultimately, we find that PSI values do not show dramatic deviation for splicing events between treatments, nor between sexes (Figure S7). The only two instances that there appears to be a binary shift is in the case of alternate acceptors (AA), specifically in the hypothalamus of the male vs female analysis, and in the hypothalamus of the control vs stress analysis. This suggests that sex is not a driving factor for PSI values of stress-response alternative splice sites.
We would like to note that the algorithm used to identify significant splice events implicitly biases our results toward a bimodal distribution as any event with a deltaPSI between − 0.1 and 0.1 was discounted. To determine how dramatically this shifted our distributions, we temporarily called all splicing events with p < = 0.05 significant. The additional 9 events this included for the male-female comparison were all found in the pituitary (3 control, 6 stress). Sixty-one events in the control-stress comparison had deltaPSI values between − 0.1 and 0.1, 58 of which were also in the pituitary (53 male, 5 female) and 3 in the female gonad. Most significant splicing events have abs(deltaPSI) values greater than 0.1 in the male-female comparison. This is perhaps not particularly surprising as splicing differences are likely more biologically effective when included or excluded at higher PSI levels due to a more dramatic shift in protein presence for a given isoform. The 53 events in the male pituitary with low abs(deltaPSI) may suggest the male pituitary responds to stress not only by splicing genic regions at higher PSI values, but also by numerous, low-level splicing events. Further analyses using additional splicing tools and datasets are needed to validate this hypothesis.
To survive an unanticipated environmental perturbation, or stressor, the body activates multiple physiological systems to mobilize needed resources, with the intent of returning the body back to a homeostatic level . The stress response is considered relatively conserved across vertebrates [73,74,75,76,77] and includes the limbic system, the central nervous system, the stress or hypothalamic-pituitary-adrenal (HPA) axis, and the thyroid or hypothalamic-pituitary-thyroid axis (HPT), all of which can directly or indirectly influence the HPG axis [77,78,79,80]. Selection should favor an optimized stress response in which an organism’s physiological and behavioral responsiveness to perturbation promotes its survival and, ultimately, reproduction. However, we know little of how stress can affect the reproductive system at the genomic level of alternative splicing, and even less so about the similarities and differences experienced by males and females. Using the rock dove model in a highly-replicated RNA-Seq study, we tested how an external environmental stressor could affect the alternative splicing response of the HPG axis in both males and females. We discovered a significant sex difference in the splicing profile of hypothalamic and pituitary tissues. When subjects underwent a restraint stress treatment, we observed a decrease in the total number of significant splicing events between sexes as compared to controls (Fig. 1). Two potential explanations for this are (1) the “convergent landscape hypothesis”: when males and females are stressed, they converge on a single splicing landscape, or (2) the “resource redirection hypothesis”: the stress response shuts down non-essential biological processes, minimizing background splicing “noise” that differentiates unstressed males and females, to redirect resources to bodily processes critical for survival.
Convergent landscape hypothesis
The decrease in alternative splicing events experienced by stressed males and females as compared to their unstressed controls may suggest that both sexes are converging on an optimal splicing landscape. Previous reports show differential gene expression (DGE) between sexes in both an unstressed control state [35, 46, 47], and in response to a stressor [7, 35, 46, 47]. Given these previous findings, coupled with ours presented here, we suggest that future study of the HPG stress response should endeavor to further characterize sex-specific and sex-convergent responses at all levels - genomic, proteomic, and methylomic. If other response mechanisms between males and females showed a similar pattern where the sexes exhibit fewer differences in a stress state, it may suggest that there was an optimal stress response phenotype being maintained at numerous levels, in addition to the genomic levels currently understood. If this pattern of convergence is not seen at other physiological levels, it is possible that the decrease in splicing differences between males and females in the restraint stress treatment may reflect a stress response unique to the splicing level that does not exist in the proteomic or methylomic levels or additional variables unassessed in this study (e.g., stage of reproductive cycle).
Resource redirection hypothesis
A second hypothesis to explain the presence of decreased alternative splicing differences between stressed males and females relative to controls is that the stress response redirects resources from less essential biological processes, like reproduction, to those more necessary for survival. Biological sex manifests in various distinct physiological phenotypes, particularly concerning those associated with reproductive processes [81,82,83,84]. However, an organism’s response to a stressor, no matter the sex, can trigger the mobilization of internal bodily resources to promote its survival [85,86,87,88]. Thus, an active sex-specific physiologically-driven alternative splicing landscape under stress may become reconciled in the HPG axis to shunt focus to other biological processes to support survival. This reduction in sexually dimorphic splicing background noise would result in fewer splicing differences between the sexes, as we observed. The remaining splicing differences may lie in loci generally implicated in the stress response, though future studies are needed to say this definitively. While our attempts to map spliced genes to Kegg pathways did not yield a discernable pattern of splicing of a particular pathway, we did find that several of the alternatively spliced genes may be involved in DNA-damage response and immune system processes, both of which are intimately linked with the stress response [89, 90]. Additional study of the splicing landscape of males and females in various conditions, environments, age groups, and reproductive stages will shed further light on splicing variation between the sexes, and if our results seen are indeed a reflection of decreased activity of peripheral biological processes during stress.
Genomic events by type
As compared to the hypothalamus and gonads, the pituitary in both sexes exhibited the most active splicing profile under stress as compared to controls. This pattern was consistent with our previous report of increased differential gene expression (DGE) in the pituitary following restraint stress , suggesting this gland may be a targeted site for mediating a stress response in the HPG axis through both gene expression and alternative splicing mechanisms. Regarding overall splicing activity, males exhibited greater responsiveness to stress than females in the hypothalamus and pituitary, with a significant increase in response of the male hypothalamus compared to that of the female. Contrarily, examination of the same samples taken from the same birds showed that females exhibited greater DGE responsiveness to stress than males in all three tissues . Thus, while examination at a gene transcription level may suggest the female reproductive axis as more responsive to the restraint stress treatment, examination at the alternative splicing level suggests males to be the more responsive of the sexes. The next steps would be to examine the protein products produced, but for now suffice it to say caution should be taken when determining which sex is more “responsive” to any stimuli, as this can vary by biological level of examination.
Further investigation of the types of splicing events that occurred following restraint stress revealed that core exon events as compared to all other event types were vastly more abundant across all tissues, regardless of sex or treatment, a phenomenon supported by findings in humans and mice . Regarding the effects of stress treatment, more core exon events were present in the female gonads as compared to the male gonads. The higher abundance of core exon events could indicate that exon inclusion/exclusion is a more efficient approach to modulate a resulting protein, thereby triggering a more specific phenotypic response to a change in environment [20, 92,93,94]. This may then signify that the stress response of the female gonad is driven more heavily by core exon alterations to a transcript, while transcripts in the male gonad are more readily modified with all types of splicing events.
The only other splicing event type that differed significantly between the stress response of males and females was the retention of introns, with more retained intron events occurring in the male hypothalamus. Retained intron events are a means by which gene expression can be modified, often by downregulation . However, these events may also alter cis and trans elements that would otherwise modify transcript stability or translational efficacy . Perhaps through differential utilization of the retained intron splicing event, the male hypothalamus could further modulate gene expression, targeting a more optimal stress response.
In an attempt to gain a better understanding of potential products being generated by splicing activity, we examined our results using Gene Ontology (GO). Our inquiry of genes undergoing sex-biased alternative splicing revealed many GO terms related to “binding.” The function of these genes is linked with how hormones bind to their receptors, and by extension, how endocrine messages are received. With the increase in abundance of spliced hypothalamic binding genes and a decrease of “binding” terms in the pituitary spliced genes, it appears that males and females may differ in how cells in the hypothalamus bind various protein products. While the function of genes in the male and female pituitary related to “binding” exhibit fewer alternative splicing events than expected by genomic prediction, the sexes showed greater difference in alternative splicing of genes associated with “cell/neuronal structure”. This may suggest that similar functions can arise from sex-specific structures.
Alternative splicing exhibited in stressed versus control treatment subject HPG tissues also reveal sexual-dimorphic functionality. Once again, many of the alternatively spliced gene products within each sex in response to stress are related to “binding”. Specifically, we observed elevated splicing of these genes in the male hypothalamus, female pituitary, and male gonad. This pattern seems to suggest that the female pituitary is responding to restraint stress by altering how molecules/drugs are bound, while the male hypothalamus is more active in responding to stress via splicing binding genes. Thus, sexual differences in splicing landscapes in our samples may signify sex-biased hypothalamic function in response to stress, but similar pituitary functions between the sexes.
Exon size, location, and motifs
The exon size distribution of alternative splicing events was significantly smaller than that predicted by the genomic distribution, agreeing with previous findings that the likelihood of a core exon being spliced may be inversely related to its size . This inverse relationship of exon size and frequency of splicing may be a reflection of the highly-specific and precise adjustments required of splicing machinery [15, 98, 99]. Inclusion or excision of a larger genetic sequence could prove more biologically challenging due to the physical distance between the two splice sites, and variable inclusion of a larger exon may result in a more dramatic alteration to the protein product, negating the effect of fine-tuning a response. Though the alternative splicing detected in both sexes exhibited smaller than expected exon sizes, size distributions in the female pituitary were smaller than those of the male pituitary. This may suggest that splicing in the female pituitary is targeted towards minute alteration, while genes in the male pituitary are modulated by movement/removal of larger sequence segments. Additional study of splicing in pituitaries of both sexes will be required to clarify the biological implications of this sex-bias.
In all comparisons, we found that the vast majority of splicing occurs in the protein coding region, or coding sequence (CDS) of a gene. When we examine the frequency of splicing location in relation to the genomic expectation (assuming an even distribution of splicing across all genic regions), we see that variation from expectation actually lies in the untranslated regions. Analyses comparing males and females, within each treatment group, revealed increased splicing of 5′-UTR genes in both the hypothalamus and pituitary, and a > 6% decrease in splicing of the 3′-UTR genes in the pituitary. Splicing differences identified within the HPG axis of each sex in stressed birds as compared to controls displayed lower 3’UTR events in all tissues, and lower 5’UTR events in the pituitary. The increased abundance of 5’UTR alternative splicing events in males as compared to female subjects suggests that differences in the 5’UTR region of transcripts exist between the sexes (regardless of control/stress state). Previous research has found that the 5’UTR can mediate translation of proteins, regulate expression, and even regulate some splicing events [100,101,102]. Thus, differential splicing of the 5’UTR between males and females may be a sex-specific mechanism for modifying protein production. The one instance in which we see a decrease in splicing of 5’UTR events is the pituitary, where each sex exhibits fewer splicing events in response to stress. Males and females may differ in the 5’UTR, but both sexes seem to minimize splicing of this region in the pituitary in response to restraint stress. The decrease in splicing events in the 3’UTR could signify that splicing events in the HPG axis do not favor the 3’UTR in response to restraint stress. 3′ ends of transcripts may be protected from alternative splicing, as evidenced by this decrease observed in both sexes and treatments. Degradation of mRNA transcripts can increase due to environmental stress , and it may be that splicing events in the 3’UTR would result in undesired degradation of transcripts, or perhaps no degradation at all.
We did not identify any differences in the frequency of protein motif splicing between the sexes or treatments. As motifs combinations and conformations heavily impact protein function, we anticipated identifying specific motifs wherein splicing more readily occurred. Our query did not reveal this, suggesting that splicing is not targeting protein motifs within the HPG axis in either sex or treatment. Instead, the genic region (CDS, 5’UTR, 3’UTR) and gene-type appear to be larger drivers of the location that splicing events occur.
In both sexes and treatment groups, we identified more alternative splicing activity in the pituitary as compared to the hypothalamus and gonads. Our previous research also revealed the pituitary to have an increased response to restraint stress at the level of gene expression as compared to the hypothalamus or gonads , suggesting a heightened stress response in this part of the HPG axis. Of genomic splicing events identified, core exon events are the most abundant form of splicing in all tissues/sexes/treatments, and splicing more frequently occurs in CDS regions. We found less splicing in the 3’UTR and more splicing in the 5’UTR than expected in both sex and treatment comparisons, suggesting that the 3′ end of transcripts is more biologically constrained than the 5′ end in our samples. Another splicing constraint evidenced by our findings, is that of an inverse relationship of exon size and splicing frequency. The overall reduction in splicing differences between the sexes when animals experienced a restraint stress treatment may point to a conserved splicing response to stress within the HPG axis. However, despite this reduction in sex-splicing differences, the male hypothalamus and female ovary experienced increased splicing activity in the face of stress as compared to the female hypothalamus and male testes, respectively. Sex differences in alternative splicing within the HPG axis, as well as in previously reported gene expression  support sex-specific mechanisms for the stress response of the reproductive axis. While females experience increased differential gene expression in their HPG axis as compared to males in response to restraint stress , we found that males experience increased alternative splicing. Our examination of the vertebrate stress response at multiple levels of biological organization offers a more complete picture of its mechanistic underpinnings between the sexes at an unprecedented proximate level. These data inspire further integrative levels of investigation to inform and potentially revolutionize evolutionary theory and the study of stress-induced reproductive dysfunction.
Organisms for this study were obtained and sampled at the University of California (UC), Davis where the Calisi lab maintains a breeding colony of Rock Doves. Doves were housed in semi-enclosed aviaries (5′ × 4′ × 7′), with 8 sexually reproductive adult pairs per aviary. Food and water were provided ad libitum. Birds were exposed to natural light, which was supplemented with artificial light on a 14 h of light:10 h of darkness cycle. We sampled sexually mature males and females that were without eggs or chicks so as to control for physiological changes that occur to facilitate parental care behaviors. Birds were sampled between 0900—1100 (Pacific Standard Time) to also control for potential circadian rhythm confounds. All handling and sampling procedures followed approved animal care and handling protocols (UC Davis Institutional Animal Care and Use Committee #18895). Samples from the control group were taken within 5 min of entering their cages. Birds in the restraint stress treatment group were captured in < 1 min upon entering their aviary and immediately and individually restrained in cloth bags for 30 min prior to sampling, replicating methods described in .
Prior to sampling, birds were anesthetized using isoflurane – a standard methodology amongst avian scientists that is approved by the Institutional Animal Care and Use Committee. Briefly, birds were directly exposed to cotton balls moistened with isoflurane in a 50 ml falcon tube until unconscious, at which point they were immediately decapitated. Trunk blood was collected for the detection of concentrations of the adrenal hormone, corticosterone; high concentrations confirmed the effectiveness of the stress treatment protocol . All tissues (brains, pituitaries, and gonads) were flash frozen on dry ice and transferred to a -80 °C freezer for storage until additional processing, a process described in [7, 35, 46, 47]. Briefly, a cryostat (Leica CM model 1860) was used to section the brains coronally at 100 μm to enable optimal visualization and biopsy of the hypothalamus, which included the adjoined lateral septum. Hypothalamic sections, pituitaries, and gonads were preserved in RNALater (Invitrogen, Thermo Fisher Scientific) and shipped on dry ice from UC Davis to the University of New Hampshire for cDNA preparation and sequencing. We sequenced tissue from whole homogenized hypothalami, pituitaries, testes and (homogenized tissue from the oviduct and ovarian follicles). In total, we processed and sequenced hypothalami, pituitary glands, and gonads (testes/ovaries) from 48 birds (12 males and 12 females, each, per treatment —control and stress), resulting in 144 cDNA libraries.
cDNA library preparation and sequencing
Library preparation and sequencing has been previously described in [7, 35, 46, 47]. In brief, all tissue samples were thawed on ice in an RNAse-free work environment. Total RNA was extracted with a standard Trizol extraction protocol (Thermo Fisher Scientific, Waltham, MA). Total RNA quality and quantity was assessed with gel electrophoresis and a Broad Range RNA Qubit assay (Thermo Fisher Scientific, Waltham, Massachusetts). Illumina sequence libraries were prepared using the NEB Next Ultra Directional RNA Library Prep Kit (New England Biolabs, Ipswich, MA). Quality and quantity of the cDNA libraries were validated with a Tapestation 2200 Instrument (Agilent, Santa Clara, California). Libraries were diluted to 5 nM and pooled. We sent multiplexed library pools to the New York Genome Center for 125 base pair paired-end sequencing on a HiSeq 2500 platform. Sixteen samples were re-sequenced at Novogene Sequencing Center to increase coverage of those samples with the lowest read abundance.
Pigeon genome annotation
Many of the currently-available methods for identification of splicing events are heavily dependent upon genome assembly completeness and annotation quality. The existing annotated assembly for C. livia (GCF_000337935.1) did not provide the level of contiguity required for accurate identification of splicing events. To address this, we used an unannotated chromosome-level assembly for the rock dove (GCA_001887795.1). We annotated this chromosome-level assembly with Maker version 2.3.1 . A custom script (available here: https://github.com/AndrewLangvt/Dissertation/blob/master/cliv2_gff2gtf.py) was used to generate the required gene transfer format (GTF) file from the resultant gene feature format (GFF).
Sequence QC, read processing, and splicing event identification
Sequence data were downloaded to the Pittsburg Supercomputing Center network “Bridges” and read quality confirmed with FastQC v0.11.5 (Andrews, 2011). Read data from all 144 samples was first error corrected with Rcorrector v3 , then adapter sequences and corrected reads with Phred< 2 were removed from sequence datasets.
To accurately identify splicing events in our study system, we used Whippet, an algorithm that rapidly identifies splicing events from RNAseq data . (Code available here: https://github.com/AndrewLangvt/Scripts/tree/master/splicing_analysis/whippet).
Whippet uses annotated transcript features to generate splice graphs and transcript-level mapping data to calculate Percent Spliced In (PSI) values for all paths through the directed graph . Whippet v0.10.4 was used to quantify read data and identify differential splicing events between treatments. Events with probability > 0.95 and abs(deltaPSI) ≥0.1 were deemed significant and considered to be true splicing events.
Splicing event comparisons
To understand the sex-specific response to restraint stress, we asked (1) how similar are males and females at the exonic level (i.e. assessing splicing events between sexes, within treatments), and (2) how do males and females differ in their response to stress (i.e., assessing splicing events within the sexes, between treatments). We did not include analyses between different tissues (e.g., male gonads to female gonads, or male pituitary to male hypothalamus) as numerous other studies have identified tissue-specific splicing events; these comparisons would not aid in discerning unique splicing patterns between sexes, or between treatments [107,108,109]. Whippet identifies splicing events by type, including Core Exon (CE), Alternate Acceptor (AA), Alternate Donor (AD), Retained Intron (RI), Tandem Alternative Polyadenylation Site (TE), Tandem Transcription Start Site (TS), Alternative First Exon (AF), and Alternative Last Exon (AL). Of all the event types detected by Whippet, core exon, alternate acceptor, alternate donor, and retained intron encompassed 95% of our identified splicing events. For simplicity, we did not represent the TE, TS, AF, or AL events in our figure depicting splicing by type even though these events are included in all of the analyses. Whippet quantifies splicing abundance to generate percent spliced in (PSI) values for each splice event by sample. Next, PSI values across all sites for samples of each treatment group are compared via the “delta” script provided with this tool. Our high-replicate study design enhances the power of this tool in accurately calling AS events. To increase stringency of calling AS events, Whippet requires perfect 18 nt matches over exon-exon junctions, thus limiting multi-mapping events. To further ensure we were correctly identifying splicing events, we filtered out events of low probability and low PSI. From the total number of events identified by Whippet, only those with a probability > 0.95 and an absolute value of PSI > 0.1 were considered real events. This method of filtering allowed us to retain only true positives in our results.
Gene ontology parent term analysis
We assessed lists of genes exhibiting AS for gene ontology (GO) term enrichment, though no statistically significant term enrichments were found. As such, we characterized our lists of spliced genes by identifying abstracted (or parent) GO terms of spliced genes to identify patterns of enrichment or depletion of overarching gene functions. Lists of Entrez Gene IDs were analyzed through the “Gene List Analysis” tool at http://www.pantherdb.org to link GO terms with each ID. These terms were further analyzed with go_abstract.py (https://github.com/twestbrookunh), a python script that abstracts a provided ontology term to its “parent” term from the Gene Ontology Consortium Database [110,111,112]. We generated a matrix of expected abundances for each tissue using this same analysis on the entire C. livia genome. GO terms with counts deemed significantly different from expected counts via Chi-Squared analyses were visualized in R. We calculated differences of counts from expected values and did not include parent ontology terms of very low abundance in the genome (< 2%) in our visualization.
Exon size, location, and motif analysis
We used a custom Python script, available at https://github.com/AndrewLangvt/Scripts/blob/master/splicing_analysis/exon_extract.py to extract desired exons from the genome and translate them into amino acid sequences based upon the phase and strand delineated in the GFF while simultaneously determining the length of each exon and the transcript region from which it was spliced (5′-UTR, CDS- coding sequence, or 3′-UTR). Distributions of exon lengths were compared to the genomic distribution to determine statistical significance using a Wilcoxon test. Counts of splicing event by region (5′-UTR, CDS, or 3′-UTR) were normalized by the total number of possible splicing locations within each region, thus negating any differential in abundance of potential splice sites by region. Variation from expected abundance of splicing event by region was then normalized by expected values across tissues. Expected counts of splicing events per region were determined from the genomic proportions and difference from expected was normalized by expected values across tissues. We further assessed the functional nature of alternatively spliced exons by investigating these sequences for protein motifs. To accomplish this, we used MOTIF, a GenomeNet tool (https://www.genome.jp/tools/motif/) which identifies protein motifs within protein sequences by searching the PFAM database  and NCBI databases COG (PMID: 10592175), SMART (PMID: 10592234), and TIGRFAM (PMID: 11125044).
Availability of data and materials
Read data corresponding to the control birds are available using the European Nucleotide Archive project ID PRJEB16136; read data corresponding to the stressed birds are available at PRJEB21082. All code for analyses in this manuscript can be found at https://github.com/AndrewLangvt/Scripts/tree/master/splicing_analysis/.
Alternative First Exon
Alternative Last Exon
complementary deoxyribonucleid acid
Differential gene expression
Gene feature format
Gene transfer format
Percent Spliced In
Tandem Alternative Polyadenylation Site
Tandem Transcription Start Site
University of California
Holway DA, Suarez AV. Animal behavior: an essential component of invasion biology. Trends Ecol Evol. 1999;14:328–30.
Sutherland WJ. The importance of behavioural studies in conservation biology. Anim Behav. 1998;56:801–9.
Vermeij GJ. When biotas meet: understanding biotic interchange. Science. 1991;253:1099–104.
Wong BBM, Candolin U. Behavioral responses to changing environments. Behav Ecol. 2015;26:665–73.
Sapolsky RM, Romero LM, Munck AU. How do glucocorticoids influence stress responses? Integrating permissive, suppressive, stimulatory, and preparative actions. Endocr Rev. 2000;21:55–89.
Zhao C, Wang P, Si T, Hsu C-C, Wang L, Zayed O, Yu Z, Zhu Y, Dong J, Tao WA, Zhu J-K. MAP kinase cascades regulate the cold response by modulating ICE1 protein stability. Dev Cell. 2017;43:618–629.e5.
Calisi RM, Austin SH, Lang AS, MacManes MD. Sex-biased transcriptomic response of the reproductive axis to stress. Horm Behav. 2018;100:56–68.
Ohama N, Sato H, Shinozaki K, Yamaguchi-Shinozaki K. Transcriptional regulatory network of plant heat stress response. Trends Plant Sci. 2017;22:53–65.
Smith CW, Patton JG, Nadal-Ginard B. Alternative splicing in the control of gene expression. Annu Rev Genet. 1989;23:527–77.
Tapial J, Ha KCH, Sterne-Weiler T, Gohr A, Braunschweig U, Hermoso-Pulido A, Quesnel-Vallières M, Permanyer J, Sodaei R, Marquez Y, Cozzuto L, Wang X, Gómez-Velázquez M, Rayon T, Manzanares M, Ponomarenko J, Blencowe BJ, Irimia M. An atlas of alternative splicing profiles and functional associations reveals new regulatory programs and genes that simultaneously express multiple major isoforms. Genome Res. 2017;27:1759–68.
Izquierdo J-M, Valcárcel J. A simple principle to explain the evolution of pre-mRNA splicing. Genes Dev. 2006;20:1679–84.
Gudlaugsdottir S, Boswell DR, Wood GR, Ma J. Exon size distribution and the origin of introns. Genetica. 2007;131:299–306.
Krawczak M, Thomas NST, Hundrieser B, Mort M, Wittig M, Hampe J, Cooper DN. Single base-pair substitutions in exon-intron junctions of human genes: nature, distribution, and consequences for mRNA splicing. Hum Mutat. 2007;28:150–8.
Sakharkar MK, Chow VTK, Kangueane P. Distributions of exons and introns in the human genome. In Silico Biol. 2004;4:387–93.
Birzele F, Csaba G, Zimmer R. Alternative splicing and protein structure evolution. Nucleic Acids Res. 2008;36:550–8.
Bush SJ, Chen L, Tovar-Corona JM, Urrutia AO. Alternative splicing and the evolution of phenotypic novelty. Philos Trans R Soc Lond B Biol Sci. 2017;372. https://doi.org/10.1098/rstb.2015.0474.
Kim E, Magen A, Ast G. Different levels of alternative splicing among eukaryotes. Nucleic Acids Res. 2007;35:125–31.
Gueroussov S, Gonatopoulos-Pournatzis T, Irimia M, Raj B, Lin Z-Y, Gingras A-C, Blencowe BJ. An alternative splicing event amplifies evolutionary differences between vertebrates. Science. 2015;349:868–73.
Tian B, Manley JL. Alternative polyadenylation of mRNA precursors. Nat Rev Mol Cell Biol. 2017;18:18–30.
Wang Y, Liu J, Huang BO, Xu Y-M, Li J, Huang L-F, Lin J, Zhang J, Min Q-H, Yang W-M, Wang X-Z. Mechanism of alternative splicing and its regulation. Biomed Rep. 2015;3:152–8.
Wang Z, Burge CB. Splicing regulation: from a parts list of regulatory elements to an integrated splicing code. RNA. 2008;14:802–13.
Wang Z, Xiao X, Van Nostrand E, Burge CB. General and specific functions of exonic splicing silencers in splicing control. Mol Cell. 2006;23:61–70.
Preußner M, Goldammer G, Neumann A, Haltenhof T, Rautenstrauch P, Müller-McNicoll M, Heyd F. Body temperature cycles control rhythmic alternative splicing in mammals. Mol Cell. 2017;67:433–446.e4.
Staiger D, Brown JWS. Alternative splicing at the intersection of biological timing, development, and stress responses. Plant Cell. 2013;25:3640–56.
Sun H. Deciphering alternative splicing and nonsense-mediated decay modulate expression in primary lymphoid tissues of birds infected with avian pathogenic E. coli (APEC). BMC Genet. 2017;18:21.
Laloum T, Martín G, Duque P. Alternative splicing control of abiotic stress responses. Trends Plant Sci. 2018;23:140–50.
Nijholt I, Farchi N, Kye M, Sklan EH, Shoham S, Verbeure B, Owen D, Hochner B, Spiess J, Soreq H, Blank T. Stress-induced alternative splicing of acetylcholinesterase results in enhanced fear memory and long-term potentiation. Mol Psychiatry. 2004;9:174–83.
Climente-Gonzalez H, Porta-Pardo E, Godzik A, Eyras E. The functional impact of alternative splicing in cancer. bioRxiv. 2017:076653. https://doi.org/10.1101/076653.
Chen J, Weiss WA. Alternative splicing in cancer: implications for biology and therapy. Oncogene. 2015;34:1–14.
Siegfried Z, Karni R. The role of alternative splicing in cancer drug resistance. Curr Opin Genet Dev. 2018;48:16–21.
Gibilisco L, Zhou Q, Mahajan S, Bachtrog D. Alternative splicing within and between Drosophila species, sexes, tissues, and developmental stages. PLoS Genet. 2016;12:e1006464.
McIntyre LM, Bono LM, Genissel A, Westerman R, Junk D, Telonis-Scott M, Harshman L, Wayne ML, Kopp A, Nuzhdin SV. Sex-specific expression of alternative transcripts in Drosophila. Genome Biol. 2006;7:R79.
Telonis-Scott M, Kopp A, Wayne ML, Nuzhdin SV, McIntyre LM. Sex-specific splicing in Drosophila: widespread occurrence, tissue specificity and evolutionary conservation. Genetics. 2009;181:421–34.
Trabzuni D, Ramasamy A, Imran S, Walker R, Smith C, Weale ME, Hardy J, Ryten M, North American Brain Expression Consortium. Widespread sex differences in gene expression and splicing in the adult human brain. Nat Commun. 2013;4:2771.
MacManes MD, Austin SH, Lang AS, Booth A, Farrar V, Calisi RM. Widespread patterns of sexually dimorphic gene expression in an avian hypothalamic-pituitary-gonadal (HPG) axis. Sci Rep. 2017;7:45125.
Nozaki M. Hypothalamic-pituitary-gonadal endocrine system in the hagfish. Front Endocrinol. 2013;4:200.
Sower SA, Freamat M, Kavanaugh SI. The origins of the vertebrate hypothalamic-pituitary-gonadal (HPG) and hypothalamic-pituitary-thyroid (HPT) endocrine systems: new insights from lampreys. Gen Comp Endocrinol. 2009;161:20–9.
Geraghty AC, Kaufer D. Glucocorticoid Regulation of Reproduction. Adv Exp Med Biol. 2015;872:253–78.
Johnson JE, Kalmar GB, Sohal PS, Walkey CJ, Yamashita S, Cornell RB. Comparison of the lipid regulation of yeast and rat CTP: phosphocholine cytidylyltransferase expressed in COS cells. Biochem J. 1992;285(Pt 3):815–20.
Toufexis D, Rivarola MA, Lara H, Viau V. Stress and the reproductive axis. J Neuroendocrinol. 2014;26:573–86.
Ball GF, Silver R. Timing of incubation bouts by ring doves (Streptopelia risoria). J Comp Psychol. 1983;97:213–25.
Buntin JD, Ruzycki E, Witebsky J. Prolactin receptors in dove brain: autoradiographic analysis of binding characteristics in discrete brain regions and accessibility to blood-borne prolactin. Neuroendocrinology. 1993;57:738–50.
Lehrman DS. The physiological basis of parental feeding behavior in the ring dove (Streptopelia Risoria). Behaviour. 1955;7:241–85.
Darwin C. On the origin of species by means of natural selection, or the preservation of favoured races in the struggle for life. British Foreign Medico Chirurgical Rev. 1860;25:367–404.
Domyan ET, Kronenberg Z, Infante CR, Vickrey AI, Stringham SA, Bruders R, Guernsey MW, Park S, Payne J, Beckstead RB, Kardon G, Menke DB, Yandell M, Shapiro MD. Molecular shifts in limb identity underlie development of feathered feet in two domestic avian species. eLife. 2016;5:e12115.
Gillespie MJ, Crowley TM, Haring VR, Wilson SL, Harper JA, Payne JS, Green D, Monaghan P, Stanley D, Donald JA, Nicholas KR, Moore RJ. Transcriptome analysis of pigeon milk production - role of cornification and triglyceride synthesis genes. BMC Genomics. 2013;14:169.
Shapiro MD, Kronenberg Z, Li C, Domyan ET, Pan H, Campbell M, Tan H, Huff CD, Hu H, Vickrey AI, Nielsen SCA, Stringham SA, Hu H, Willerslev E, Gilbert MTP, Yandell M, Zhang G, Wang J. Genomic diversity and evolution of the head crest in the rock pigeon. Science. 2013;339:1063–7.
Chandran UR, DeFranco DB. Regulation of gonadotropin-releasing hormone gene transcription. Behav Brain Res. 1999;105:29–36.
Cheng CK, Yeung CM, Hoo RLC, Chow BKC, Leung PCK. Oct-1 is involved in the transcriptional repression of the gonadotropin-releasing hormone receptor gene. Endocrinology. 2002;143:4693–701.
Campbell BK, Dobson H, Scaramuzzi RJ. Ovarian function in ewes made hypogonadal with GnRH antagonist and stimulated with FSH in the presence or absence of low amplitude LH pulses. J Endocrinol. 1998;156:213–22.
Pohl CR, Richardson DW, Hutchison JS, Germak JA, Knobil E. Hypophysiotropic signal frequency and the functioning of the pituitary-ovarian system in the rhesus monkey. Endocrinology. 1983;112:2076–80.
Takahashi A, Kanda S, Abe T, Oka Y. Evolution of the hypothalamic-pituitary-gonadal Axis regulation in vertebrates revealed by knockout Medaka. Endocrinology. 2016;157:3994–4002.
Jin J-M, Yang W-X. Molecular regulation of hypothalamus-pituitary-gonads axis in males. Gene. 2014;551:15–25.
Rosen H, Jameel ML, Barkan AL. Dexamethasone suppresses gonadotropin-releasing hormone (GnRH) secretion and has direct pituitary effects in male rats: differential regulation of GnRH receptor and gonadotropin responses to GnRH. Endocrinology. 1988;122:2873–80.
Schneider F, Tomek W, Gründker C. Gonadotropin-releasing hormone (GnRH) and its natural analogues: a review. Theriogenology. 2006;66:691–709.
Nellore A, Paziana K, Ma C, Tsygankova OM, Wang Y, Puttaswamy K, Iqbal AU, Franks SR, Lv Y, Troxel AB, Feldman MD, Meinkoth JL, Brose MS. Loss of Rap1GAP in papillary thyroid cancer. J Clin Endocrinol Metab. 2009;94:1026–32.
Katagiri K, Hattori M, Minato N, Kinashi T. Rap1 functions as a key regulator of T-cell and antigen-presenting cell interactions and modulates T-cell responses. Mol Cell Biol. 2002;22:1001–15.
González-Tajuelo R, Silván J, Pérez-Frías A, de la Fuente-Fernández M, Tejedor R, Espartero-Santos M, Vicente-Rabaneda E, Juarranz Á, Muñoz-Calleja C, Castañeda S, Gamallo C, Urzainqui A. P-Selectin preserves immune tolerance in mice and is reduced in human cutaneous lupus. Sci Rep. 2017;7:41841.
Matsumoto M, Miyasaka M, Hirata T. P-selectin glycoprotein ligand-1 negatively regulates T-cell immune responses. J Immunol. 2009;183:7204–11.
Ariazi EA, Cunliffe HE, Lewis-Wambi JS, Slifker MJ, Willis AL, Ramos P, Tapia C, Kim HR, Yerrum S, Sharma CGN, Nicolas E, Balagurunathan Y, Ross EA, Jordan VC. Estrogen induces apoptosis in estrogen deprivation-resistant breast cancer through stress responses as identified by global gene expression across time. Proc Natl Acad Sci U S A. 2011;108:18879–86.
Acevedo-Rodriguez A, Kauffman AS, Cherrington BD, Borges CS, Roepke TA, Laconi M. Emerging insights into hypothalamic-pituitary-gonadal axis regulation and interaction with stress signalling. J Neuroendocrinol. 2018;30:e12590.
Handa RJ, Burgess LH, Kerr JE, O’Keefe JA. Gonadal steroid hormone receptors and sex differences in the hypothalamo-pituitary-adrenal axis. Horm Behav. 1994;28:464–76.
Chinery R, Brockman JA, Dransfield DT, Coffey RJ. Antioxidant-induced nuclear translocation of CCAAT/enhancer-binding protein β: a critical role for protein kinase a-mediated phosphorylation of Ser299. J Biol Chem. 1997;272:30356–61.
Pless O, Kowenz-Leutz E, Knoblich M, Lausen J, Beyermann M, Walsh MJ, Leutz A. G9a-mediated lysine methylation alters the function of CCAAT/enhancer-binding protein-beta. J Biol Chem. 2008;283:26357–63.
Miller AL, Bement WM. Regulation of cytokinesis by rho GTPase flux. Nat Cell Biol. 2009;11:71–7.
Moon SY, Zheng Y. Rho GTPase-activating proteins in cell regulation. Trends Cell Biol. 2003;13:13–22.
Godoy J, Nishimura M, Webster NJG. Gonadotropin-releasing hormone induces miR-132 and miR-212 to regulate cellular morphology and migration in immortalized LbetaT2 pituitary gonadotrope cells. Mol Endocrinol. 2011;25:810–20.
Naor Z, Benard O, Seger R. Activation of MAPK cascades by G-protein-coupled receptors: the case of gonadotropin-releasing hormone receptor. Trends Endocrinol Metab. 2000;11:91–9.
Gerhardt C, Lier JM, Burmühl S, Struchtrup A, Deutschmann K, Vetter M, Leu T, Reeg S, Grune T, Rüther U. The transition zone protein Rpgrip1l regulates proteasomal activity at the primary cilium. J Cell Biol. 2015;210:115–33.
Tokue S-I, Sasaki M, Nakahata N. Thromboxane A2-induced signal transduction is negatively regulated by KIAA1005 that directly interacts with thromboxane A2 receptor. Prostaglandins Other Lipid Mediat. 2009;89:8–15.
Struchtrup A, Wiegering A, Stork B, Rüther U, Gerhardt C. The ciliary protein RPGRIP1L governs autophagy independently of its proteasome-regulating function at the ciliary base in mouse embryonic fibroblasts. Autophagy. 2018;14:567–83.
Michael Romero L, Wingfield JC. Tempests, poxes, predators, and people: stress in wild animals and how they cope. United Kingdom: Oxford University Press; 2015.
Clarke AS, Wittwer DJ, Abbott DH, Schneider ML. Long-term effects of prenatal stress on HPA axis activity in juvenile rhesus monkeys. Dev Psychobiol. 1994;27:257–69.
Gutierrez-Triana JA, Herget U, Lichtner P, Castillo-Ramírez LA, Ryu S. A vertebrate-conserved cis-regulatory module for targeted expression in the main hypothalamic regulatory region for the stress response. BMC Dev Biol. 2014;14:41.
Leonard BE. The HPA and immune axes in stress: the involvement of the serotonergic system. Eur Psychiatry. 2005;20:S302–6.
Taniuchi S, Miyake M, Tsugawa K, Oyadomari M, Oyadomari S. Integrated stress response of vertebrates is regulated by four eIF2α kinases. Sci Rep. 2016;6:32886.
Tsigos C, Chrousos GP. Hypothalamic–pituitary–adrenal axis, neuroendocrine factors and stress. J Psychosom Res. 2002;53:865–71.
Helmreich DL, Parfitt DB, Lu X-Y, Akil H, Watson SJ. Relation between the hypothalamic-pituitary-thyroid (HPT) axis and the hypothalamic-pituitary-adrenal (HPA) axis during repeated stress. Neuroendocrinology. 2005;81:183–92.
Herman JP, Ostrander MM, Mueller NK, Figueiredo H. Limbic system mechanisms of stress regulation: hypothalamo-pituitary-adrenocortical axis. Prog Neuro Psychopharmacol Biol Psychiatry. 2005;29:1201–13.
Romero-Ramírez L, Nieto-Sampedro M, Barreda-Manso MA. Integrated stress response as a therapeutic target for CNS injuries. Biomed Res Int. 2017;2017:6953156.
Arnold AP, Itoh Y. Factors causing sex differences in birds. Avian Biol Res. 2011;4. https://doi.org/10.3184/175815511X13070045977959.
Bundy JL, Vied C, Nowakowski RS. Sex differences in the molecular signature of the developing mouse hippocampus. BMC Genomics. 2017;18:237.
Lee SK. Sex as an important biological variable in biomedical research. BMB Rep. 2018;51:167–73.
Ngun TC, Ghahramani N, Sánchez FJ, Bocklandt S, Vilain E. The genetics of sex differences in brain and behavior. Front Neuroendocrinol. 2011;32:227–46.
Goldstein DS. Adrenal responses to stress. Cell Mol Neurobiol. 2010;30:1433–40.
Nicolaides NC, Kyratzi E, Lamprokostopoulou A, Chrousos GP, Charmandari E. Stress, the stress system and the role of glucocorticoids. Neuroimmunomodulation. 2015;22:6–19.
Smulders TV. The avian hippocampal formation and the stress response. Brain Behav Evol. 2017;90:81–91.
Wendelaar Bonga SE. The stress response in fish. Physiol Rev. 1997;77:591–625.
Dhabhar FS. Effects of stress on immune function: the good, the bad, and the beautiful. Immunol Res. 2014;58:193–210.
Hara MR, Kovacs JJ, Whalen EJ, Rajagopal S, Strachan RT, Grant W, Towers AJ, Williams B, Lam CM, Xiao K, Shenoy SK, Gregory SG, Ahn S, Duckett DR, Lefkowitz RJ. A stress response pathway regulates DNA damage through β2-adrenoreceptors and β-arrestin-1. Nature. 2011;477:349–53.
Koralewski TE, Krutovsky KV. Evolution of exon-intron structure and alternative splicing. PLoS One. 2011;6:e18055.
Chen M, Manley JL. Mechanisms of alternative splicing regulation: insights from molecular and genomics approaches. Nat Rev Mol Cell Biol. 2009;10:741–54.
Rodrigues R, Grosso AR, Moita L. Genome-wide analysis of alternative splicing during dendritic cell response to a bacterial challenge. PLoS One. 2013;8:e61975.
Yabas M, Elliott H, Hoyne GF. The role of alternative splicing in the control of immune homeostasis and cellular differentiation. Int J Mol Sci. 2015;17. https://doi.org/10.3390/ijms17010003.
Ge Y, Porse BT. The functional consequences of intron retention: alternative splicing coupled to NMD as a regulator of gene expression. Bioessays. 2014;36:236–43.
Jacob AG, Smith CWJ. Intron retention as a component of regulated gene expression programs. Hum Genet. 2017;136:1043–57.
Song S, Huang Q, Guo J, Li-Ling J, Chen X, Ma F. Comparative component analysis of exons with different splicing frequencies. PLoS One. 2009;4:e5387.
Black DL. Protein diversity from alternative splicing: a challenge for bioinformatics and post-genome biology. Cell. 2000;103:367–70.
Mittendorf KF, Deatherage CL, Ohi MD, Sanders CR. Tailoring of membrane proteins by alternative splicing of pre-mRNA. Biochemistry. 2012;51:5541–56.
Álvarez D, Voß B, Maass D, Wüst F, Schaub P, Beyer P, Welsch R. Carotenogenesis is regulated by 5’UTR-mediated translation of Phytoene synthase splice variants. Plant Physiol. 2016;172:2314–26.
Kramer M, Sponholz C, Slaba M, Wissuwa B, Claus RA, Menzel U, Huse K, Platzer M, Bauer M. Alternative 5′ untranslated regions are involved in expression regulation of human heme oxygenase-1. PLoS One. 2013;8:e77224.
Ohta S, Nakagawara S, Hirai S, Miyagishima K, Horiguchi G, Kodama H. The 5′ UTR intron-mediated enhancement of constitutive splicing of the tobacco microsome ω-3 fatty acid desaturase gene. Plant Biotechnol Rep. 2018;12:105–14.
Canadell D, García-Martínez J, Alepuz P, Pérez-Ortín JE, Ariño J. Impact of high pH stress on yeast gene expression: a comprehensive analysis of mRNA turnover during stress responses. Biochim Biophys Acta. 2015;1849:653–64.
Cantarel BL, Korf I, Robb SMC, Parra G, Ross E, Moore B, Holt C, Sánchez Alvarado A, Yandell M. MAKER: an easy-to-use annotation pipeline designed for emerging model organism genomes. Genome Res. 2008;18:188–96.
Song L, Florea L. Rcorrector: efficient and accurate error correction for Illumina RNA-seq reads. GigaScience. 2015;4:48.
Sterne-Weiler T, Weatheritt RJ, Best A, Ha KCH, Blencowe BJ. Whippet: an efficient method for the detection and quantification of alternative splicing reveals extensive transcriptomic complexity. bioRxiv. 2017:158519. https://doi.org/10.1101/158519.
Grosso AR, Gomes AQ, Barbosa-Morais NL, Caldeira S, Thorne NP, Grech G, von Lindern M, Carmo-Fonseca M. Tissue-specific splicing factor gene expression signatures. Nucleic Acids Res. 2008;36:4823–32.
Pineda JMB, Bradley RK. Most human introns are recognized via multiple and tissue-specific branchpoints. Genes Dev. 2018;32:577–91.
Saha A, Kim Y, Gewirtz ADH, Jo B, Gao C, McDowell IC, GTEx Consortium, Engelhardt BE, Battle A. Co-expression networks reveal the tissue-specific regulation of transcription and splicing. Genome Res. 2017;27:1843–58.
Ashburner M, Ball CA, Blake JA, Botstein D, Butler H, Cherry JM, Davis AP, Dolinski K, Dwight SS, Eppig JT, Harris MA, Hill DP, Issel-Tarver L, Kasarskis A, Lewis S, Matese JC, Richardson JE, Ringwald M, Rubin GM, Sherlock G. Gene ontology: tool for the unification of biology. The Gene Ontology Consortium. Nat Genet. 2000;25:25–9.
Mi H, Huang X, Muruganujan A, Tang H, Mills C, Kang D, Thomas PD. PANTHER version 11: expanded annotation data from gene ontology and Reactome pathways, and data analysis tool enhancements. Nucleic Acids Res. 2017;45:D183–9.
The Gene Ontology Consortium. The Gene Ontology Resource: 20 years and still GOing strong. Nucleic Acids Res. 2019;47:D330–8.
Finn RD, Coggill P, Eberhardt RY, Eddy SR, Mistry J, Mitchell AL, Potter SC, Punta M, Qureshi M, Sangrador-Vegas A, Salazar GA, Tate J, Bateman A. The Pfam protein families database: towards a more sustainable future. Nucleic Acids Res. 2016;44:D279–85.
We thank the many undergraduate and graduate students of the MacManes and Calisi Lab for their assistance on this project, including their involvement in avian husbandry, conducting the experiment, tissue collection, RNA extraction, and cDNA library preparation.
This work was funded by the National Science Foundation Division of Integrative Organismal Systems #1455960 (to RMC and MM). The funding body did not contribute to the design of the study, nor the collection, analysis, and interpretation of the data.
Ethics approval and consent to participate
This study, to include all handling and sampling procedures, was approved by the Institutional Animal Care and Use Committee (UC Davis IACUC #18895).
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Read counts and number of reads mapped for each sample. Table S2. Genome Annotation Statistics. Table S3. List of Spliced Genes. Table S4. Figure S1. Male vs. Female Biological Process GO Analysis. Figure S2. Male vs. Female Cellular Component GO Analysis. Figure S3. Male vs. Female Exon Motifs. Figure S4. Control vs. Stress Biological Process GO Analysis. Figure S5. Control vs. Stress Cellular Component GO Analysis. Figure S6. Control vs. Stress Exon Motifs. Figure S7. Distribution of Percent Spliced In (PSI) Values for the a) control vs stress analysis and also the b) male vs female analysis.
About this article
Cite this article
Lang, A.S., Austin, S.H., Harris, R.M. et al. Stress-mediated convergence of splicing landscapes in male and female rock doves. BMC Genomics 21, 251 (2020). https://doi.org/10.1186/s12864-020-6600-6