Shared and Unique Patterns Across Clades.

Within each clade we compared expression of orthologous gene groups (OGGs) between species pairs (SI Appendix, Fig. S3). We found that the mean difference in expression between species pairs was near zero across all 1,979 OGGs (SI Appendix, Fig. S3A). Intriguingly, gene expression was least variable between the two bird species (SI Appendix, Fig. S3B), which may be explained by the fact that the mating systems of these species pairs are considerably more similar than those of the other clades (Fig. 2). Specifically, both Anthus spinoletta and Prunella modularis can form pair bonds (although less frequently in P. modularis) and exhibit direct paternal care (Fig. 2 and SI Appendix, Fig. S1 and Table S1). Despite these differences in OGG expression variation across clades, clade-specific Gene Ontology (GO) analysis revealed high conservation of GO term enrichment highlighting cell communication, signaling receptor activity, and membrane proteins as consistently associated with monogamy-related expression across clades (SI Appendix, Fig. S4).

Across Vertebrates, Gene Expression Varies Concordantly Between Males of Monogamous and Nonmonogamous Species.

To assess concordance of OGG expression variation between monogamous and nonmonogamous species across vertebrate clades, we used the differential expression analysis software package DESeq2 (33). We assessed differential OGG expression multiple times with distinct evolutionary groupings from mammals to all clades where monogamous species of distinct clades were entered as interspecific replicates of monogamy (Fig. 3). We find that differences in OGG expression are generally concordant, particularly among OGGs that exhibit larger fold-differences (Fig. 3). In addition, most OGGs maintain directional concordance. For example, many of the OGGs that show increased expression (positive log₂ fold-difference) in one evolutionary group (e.g., mammals) show increased expression in other evolutionary groups (Fig. 3). These results indicate that across vertebrates monogamous species recruit a common set of OGGs despite evolutionarily independent transitions to similar mating systems. We find that, as the evolutionary frame of reference is expanded and more distantly related clades are added to the analysis, fewer OGGs retain significance. In particular, we find a large decrease in above-threshold OGGs when amniote species pairs (i.e., mammalian and avian species pairs) are included compared with only the mammalian species pairs (Fig. 3). This finding likely reflects the decreased expression variation notable in our bird species pair (SI Appendix, Fig. S3). While the observed effect of the mating system for some OGGs may be smaller at broader taxonomic scales, adding species pairs increases the statistical power. Thus, with the exception of the comparison between mammals and birds described above, the decrease in number of OGGs meeting our threshold cutoff is quite small (Fig. 3).

When we included all species pairs, we identified 123 OGGs (6.2%) associated with monogamy across vertebrates (SI Appendix, Fig. S5). We find a number of OGGs significant at one level of analysis failing to meet the significance threshold at another (Fig. 3). Many differential expression analysis approaches, including DESeq2, rely on expression variation among biological replicates to determine differential expression. This approach is limiting when biological replicates are highly variable (34) as is the case here, where species of different clades are included as interspecific replicates of monogamy. In addition to evolutionary distance, a number of biological and technical features likely generate noise in our analysis (e.g., ecological differences among species, course-grained tissue sampling, and technical variation during sequencing). To identify shared transcriptomic patterns across monogamous species and assess whether the degree of overlap in expression between monogamous species is statistically significant, we utilize the Rank-Rank Hypergeometric Overlap (RRHO) approach (34). Comparing ranked fold-differences enables discovery of OGGs that share patterns of expression among monogamous species without requiring that expression values be similar across evolutionarily distant clades. Additionally, rather than adhering to discrete thresholds to identify candidates with similar expression, RRHO identifies candidates with coordinated directional shifts in expression using a sliding threshold (34) (Fig. 4).

Using RRHO including all 1,979 OGGs, we find an enrichment of OGG overlap in the on-diagonal extremes (i.e., at high log₂ fold-differences in expression between monogamous and nonmonogamous species in both clades being compared; Fig. 4D). Most notably, we find enriched overlap of OGGs exhibiting decreased expression in monogamous species of all clades. In all comparisons, the concordant down-regulated (i.e., bottom left) quadrant of the RRHO plots has the greatest overall significance (i.e., most significant quadrant mean −log₁₀ P value), strongly indicative of a universal signature of expression among monogamous vertebrates (Fig. 4 C and D). We generally do not find enrichment at small log₂ fold-differences with a few exceptions (e.g., frogs and voles). This observation is likely due to the fact that the majority of OGGs exhibit small fold-differences (although they may be important for monogamy-related behavior). Subtle variation in expression of genes involved in clade-specific or basic cellular functions could mask potentially important similarities at this level. A consequence of this limitation of the RRHO approach is that, for OGGs with small fold-differences, a large overlap in OGGs is required to reach significant enrichment. Interestingly, clade comparisons including frogs exhibit the most enrichment in the off-diagonal (i.e., discordant) regions. Specifically, many orthologous gene groups that show decreased expression in other groups increase in the monogamous frogs, a finding consistent with the general trend of increased gene expression in the monogamous frogs (SI Appendix, Fig. S3A). This pattern is clear in all pairwise comparisons involving frogs, suggesting that the monogamy-related expression patterns in frogs are the least similar to the other clades. Comparisons including birds show the least enrichment in OGG overlap, possibly because the mating systems of the bird species in our study are relatively similar to each other, as discussed above. Importantly, however, birds still share the same overall pattern of OGG overlap in the on-diagonal extremes and show strong similarities in expression patterns in OGGs with high differential expression (e.g., Bottom Left and Top Right quadrants of birds vs. mice and birds vs. fish, Fig. 4D). Variation among clades in overall and specific patterns of OGG enrichment may reflect a difference among clades in species-specific features of monogamy. For example, male parental care is ancestral in the poison frog clade (35). Thus, biparental care associated with monogamy in our focal species evolved from male-biased care in our frog clade rather than from female-biased care as in the other clades examined here. Future analyses—including species exhibiting some but not all reproductive, parental, and agonistic phenotypes of monogamy as well as a more anatomically focused analysis—may shed light on this variation, enabling decoupling and association of gene expression patterns with more specific behavioral phenotypes.

Discovery of Candidates Underlying Monogamy.

We identified 24 candidate genes that are robustly (i.e., independently in both differential expression and RRHO analyses) associated with monogamous mating systems across vertebrates (Fig. 5, SI Appendix, Table S7 and Dataset S1). Interestingly, across lineages, genes involved in neural development, cell–cell signaling, synaptic activity (e.g., GRM6), learning and memory (e.g., DSCAM), and cognitive function are up-regulated in monogamous males. Conversely, genes involved in gene transcription (e.g., genes associated with RNA Polymerase II) and AMPA receptor regulation are down-regulated. These differences might indicate increased neural plasticity in the face of tighter transcriptional regulation in monogamous males. A more detailed discussion of these candidate genes, including a description of neurolocalization, is provided in SI Appendix, Table S7 and Dataset S1. While the majority of the candidate OGGs identified by both approaches to differential expression analysis (DESeq2 and RRHO) show increased expression in monogamous males, the greatest overlap enrichment occurs in OGGs with decreased expression in monogamous males (Fig. 4). Many of these OGGs do not survive threshold-dependent DESeq2 analysis and thus are eliminated from our list of robust candidates. Traditional candidates associated with monogamy-related behavior (e.g., AVP and OT as well as their receptors) were not identified among those consistently expressed in monogamous males, and in some species were not identified in the RNA-seq analysis (SI Appendix, Fig. S6). This finding does not indicate that these traditional candidates are unimportant in regulating monogamy-related behavior across these clades, but rather that their roles (if they are conserved across monogamous species) are not reflected in a consistent expression pattern at the coarse neuroanatomical scale that we sampled.

Phylogenetic and Ecological Correlates.

The combined discovery of monogamy-related OGG expression by differential expression analysis (Fig. 3) and enrichment of gene expression rank correlation among clades (Fig. 4) supports the hypothesis that shared gene expression patterns underlie the behavioral expression of monogamous mating systems across vertebrate clades. However, it is possible that historical and ecological features shared between clades influence similarity in neural gene expression patterns. By design, species pairs were quite similar in their ecological attributes, differing primarily in specific characteristics of their mating systems (Fig. 2 and SI Appendix, Fig. S1 and Tables S1 and S2). Thus, the selection of species pairs by clade should minimize confounding ecological factors; however, several other factors may play a role in gene expression similarity. First, divergence times between species pairs vary between ∼2.5 and 34 million years (36). Historical contingency can bias the path of evolution such that more closely related species may be more similar due to shared evolutionary history. Second, elaboration of mating systems varies among species such that comparisons between monogamous and nonmonogamous mating systems are not equivalent across the clades. For example, the bird species included in this study share a number of mating system characteristics (Fig. 2). To assess the role of evolutionary history and mating system on gene expression divergence, we compared evolutionary and mating system distances to OGG expression distance for all species pairs. Neither evolutionary distance nor mating system distance correlated with OGG expression divergence between species pairs (Fig. 6 A and B). Notably, however, we find that the birds and the frogs are the most similar in their respective transcriptomes and are also the most similar in characteristics of mating system, while at the same time being the most distantly related of all species pairs (birds: 29 Mya; frogs: 34.2 Mya; Fig. 1 and SI Appendix, Fig. S3). When expression and mating system variation attributable to phylogeny is removed (using phylogenetic independent contrasts), we find a significant relationship between neural gene expression and mating system (Fig. 6C). Even though phylogenetic relatedness and ecological attributes affect neural transcriptome similarity across species in complex ways, together these observations indicate a critical role for mating system in driving gene expression similarity in the brain.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 1.

Representative species exhibiting monogamous mating systems (including forming pair bonds, engaging in biparental care, and defending territory) were included in the analysis. To reduce clade-specific gene expression patterns and identify neural expression patterns associated with a monogamous mating system, monogamous species were paired with a closely related nonmonogamous species. Monogamous species for each clade are shown in orange, and nonmonogamous are shown in purple. Genes were grouped into OGGs using orthoMCL (44). Of the 6,125 OGGs identified across clades, the number of OGGs and genes in parentheses with expression data are shown at each tip (SI Appendix, Table S3). The number of shared OGGs and genes in parentheses (for species pairs) at each node are shown. Median divergence estimate was obtained from TimeTree (SI Appendix, Table S5) (46). Dashed branch lengths are not drawn to scale.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 2.

Mating system and ecological distances between focal monogamous and nonmonogamous species pairs. Divergence between species pairs is indicated by color intensity. Orange indicates character states that are expressed (or more elaborated) in the monogamous species. Purple characters are expressed (or more elaborated) in the nonmonogamous species. White indicates that the character state is similar in the species pair. Profiles for all species in this study are provided in SI Appendix, Fig. S1 and Tables S1 and S2.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 3.

Volcano plots indicating which of the 1,979 OGGs identified across all clades are differentially expressed at different taxonomic levels (mammals only, mammals and birds, tetrapods, all clades). Differential expression analysis was performed using DESeq2, where the monogamous (nonmonogamous) species of each clade were included as interspecific replicates of monogamy (nonmonogamy). Black circles show no differential expression at any taxonomic level. Differential expression analysis was performed on distinct evolutionary subgroups: mammals only (399 OGGs concordantly regulated), mammals and birds (155 OGGs), tetrapods (i.e., including frogs; 126 OGGs), and all four clades of vertebrates (i.e., including fishes; 123 OGGs). OGGs with a –log₁₀ P value > 1 and a log₂ fold-difference less than −1 (blue) or greater than 1 (red), respectively, are highlighted. The darker each circle, the more concordant across clades is the expression of the OGG that it represents. As more lineages are added to the analysis, more OGGs that are significant in one analysis fall below the significance threshold in another; however, adding species pairs increases the statistical power due to the increased number of interspecific replicates; thus, with the exception of evolutionary-subgroup mammals only versus mammals and birds, the decrease in number of OGGs meeting our threshold cutoff is small.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 4.

RRHO of monogamy-related log₂ fold-differences in gene expression for the 1,979 OGGs identified across all clades. Ranked log₂ fold-differences in monogamous vs. nonmonogamous mRNA levels are binned into 44 sets of 45 OGGs from the most down-regulated to the most up-regulated in the monogamous species of each clade. OGG set overlap is compared in four quadrants defined by the transition between down- and up-regulation in each clade (A, dashed lines). The color of each pixel of the matrix (A, red square) indicates the enrichment in OGG set overlap at and above that differential expression threshold (B) and is expressed as the negative log₁₀ of the Benjamini–Yekutieli-corrected P value. Significance of the enrichment is indicated by the pixel color with warm colors indicating increased enrichment. For each pairwise comparison of clades, the strength of OGG set overlap is summarized as the most significant quadrant mean negative log₁₀ of the BY-corrected (C). Mean, median, and maximum P values for each quadrant are provided in SI Appendix, Table S6. Arrows next to the silhouettes indicate the directionality in lineage A (first) and lineage B (second) (C). RRHO analyses are shown for each pairwise clade comparison (D). Negative log₁₀ of the BY-corrected P value color scale varies across plots. Dashed lines indicate the position of the switch point from down- to up-regulation in the monogamous species of each clade. Arrows on the color scale indicate the color at P value = 0.05.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 5.

Discovery of monogamy candidate genes (OGGs). To ensure a rigorous and conservative approach toward identifying OGGs robustly associated with monogamy across clades, any given gene/OGG had to fulfill two selection criteria. First, an OGG had to show a ±1 log₂ fold expression difference between the monogamous and nonmonogamous species in at least four clades. Second, an OGG had to be among the most up- or the most down-regulated in 6 of the 10 RRHO analyses (as shown in Fig. 4); this again required an OGG to be concordantly expressed in at least four of the five lineages. Which clade (if any) is nonconcordant can vary for each OGG. Relative expression levels of these candidates are illustrated in the heatmap. Note that no clade is more often discordant than any other (SI Appendix, Table S7).

• Download figure
• Open in new tab
• Download powerpoint

Fig. 6.

Evolutionary divergence (A) and divergence in mating system characteristics (B) within clades does not correlate with gene expression distance (Spearman ρ = −0.8, P = 0.95; ρ = 0.7, P = 0.12, respectively). However, when expression variation due to phylogeny is removed (using phylogenetically independent contrasts) (C), we find a significant relationship between mating system and neural gene expression across clades (linear regression r² = 0.6, t = 3.8, P = 0.005). Mating system score was calculated as the sum of the mating system characteristics where higher values indicate more elaborated monogamy (i.e., males form pair bonds, provide direct and indirect parental care, exhibit high levels of territoriality, and are less sexually dimorphic; Fig. 2 and SI Appendix, Fig. S1 and Tables S1 and S2). A principal component analysis was performed on the differentially expressed (±1 log₂ fold difference in at least one clade) and variable (Upper quartile of variance across all species) OGGs. This subset of 401 OGGs was also used for OGG expression divergence in A and B. See SI Appendix.

Conclusions.

Using a comparative transcriptomics approach, we asked whether independent transitions to a monogamous mating system across four major clades of vertebrates are associated with shared neural gene expression patterns. A shared mechanistic basis of social behavior across distantly related clades has been documented at the level of neural circuitry where brain-region-specific expression of neurochemical genes is remarkably conserved in the Social Decision Making Network of the vertebrate fore- and midbrain (37). Further neural gene expression comparisons of aggressive behavior in bees, stickleback, and mice (4) provide support for the hypothesis that broad similarities in social behavior across species may result from independent recruitment of ancient gene modules. However, inference from this analysis was limited by the complexity of homology inference among brains in bees and brain regions sampled in mice and stickleback. We expand on this previous comparative neural gene expression study of social behavior by comparing equivalent brain dissections and providing a phylogenetically informed assessment of transcriptome similarity in monogamous species of distinct and distantly related clades. Furthermore, by carefully scoring mating system characteristics and ecological attributes and by accounting for phylogenetic nonindependence, we test alternative hypotheses underlying the origins of transcriptomic similarities that originated in independent evolutionary transitions to monogamy across vertebrates. While our sampling of whole fore- and midbrains provides a coarse-grained view of the neuromolecular correlates of mating system, this approach circumvents potential problems with the sometimes tenuous homology inference of brain regions, ensuring consistent sampling across clades.

Our analysis reveals that monogamous species of distantly related vertebrate clades share more gene expression similarities than would be expected by chance, particularly in OGGs that are highly differentially expressed in monogamous species. Shared OGG expression across these clades provides evidence of deeply homologous mechanisms underlying the suite of reproductive, parental, and agonistic behavioral phenotypes that accompany monogamous mating systems. Our discovery adds to the growing body of research highlighting shared mechanistic bases of independently evolved morphological and physiological phenotypes (12). The present work illustrates that the independent evolution of complex behavioral phenotypes, like monogamy, is better conceptualized as the product of both parallel and convergent processes where many components of the underlying mechanism arise via parallel recruitment of deeply shared gene networks and some mechanisms evolve de novo (38⇓–40). We argue that a conceptual shift is critical for achieving progress toward testing mechanistic hypotheses of the evolution of complex social behavior using comparative ‘omic approaches (as has been suggested for morphological traits; cf. ref. 41). Because the preponderance of evidence supports ancient, deeply homologous molecular and developmental pathways in evolutionary diversification, similar phenotypes should be expected to have mechanistic similarities regardless of evolutionary origins. Under this model, independent evolutionary transitions to similar behavioral phenotypes are expected to share similar mechanisms as a consequence of (i) integration of new traits with shared, ancestral physiological processes and their associated developmental and genetic pathways; and (ii) similar challenges to basic physiology that organisms experience in the evolution of these similar phenotypes. Whereas previous accounts of transcriptomic similarities among evolutionarily distinct complex behavior have focused on gene expression similarities, this conceptual approach provides a context to leverage both similarities and differences in gene expression to understand the mechanistic basis of complex trait evolution. Overall, our results indicate that independent evolutionary transitions to a monogamous mating system across vertebrate clades are accompanied by similar changes in gene expression in the male fore- and midbrain. In combination with the discovery of shared transcriptomic mechanisms associated with aggression (4), learned vocalizations in birds and humans (5), and caste differentiation in hymenoptera (7), the results presented here considerably expand our understanding of how behavioral diversity evolves.