Neurobiology of pair bonding in fishes; convergence of neural mechanisms across distant vertebrate lineages

Pair bonding has independently evolved numerous times among vertebrates. The governing neural mechanisms of pair bonding have only been studied in depth in the mammalian model species, the prairie vole, Microtus ochrogaster. In this species, oxytocin (OT), arginine vasopressin (AVP), dopamine (DA), and opioid (OP) systems play key roles in signaling in the formation and maintenance of pair bonding by targeting specific social and reward-mediating brain regions. By contrast, the neural basis of pair bonding is poorly studied in other vertebrates, and especially those of early origins, limiting our understanding of the evolutionary history of pair bonding regulatory mechanisms. We compared receptor gene expression between pair bonded and solitary individuals across eight socio-functional brain regions. We found that in females, ITR and V1aR receptor expression varied in the lateral septum-like region (the Vv/Vl), while in both sexes D1R, D2R, and MOR expression varied within the mesolimbic reward system, including a striatum-like region (the Vc); mirroring sites of action in M. ochrogaster. This study provides novel insights into the neurobiology of teleost pair bonding. It also reveals high convergence in the neurochemical mechanisms governing pair bonding across actinopterygians and sarcopterygians, by repeatedly co-opting and similarly assembling deep neurochemical and neuroanatomical homologies that originated in ancestral osteithes.

Given that the neurochemical and neuroanatomical components that underpin pair bonding in M. ochrogaster originated in early vertebrates and have remained structurally and functionally conserved, it is conceivable that in at least selective cases, they may have been repeatedly co-opted and similarly assembled into a converged regulatory neural network during independent transitions to pair bonding within the sub-phylum. Indeed, nonapeptide and DA systems appear to regulate pair bonding in other species that span several lineages, and appear to do so through targeting similar brain regions. In male and female marmosets, Callithrix penicillata, OT promotes while an OTR antagonist reduces affiliation during cohabitation with a prospective partner 93 . Similarly, in male and female tamarins, Saquinus Oedipus, urinary OT increases with intra-pair affiliation 94 . In zebra finches, T. guttata, both i.c.v. and peripheral OTR antagonist administration impairs pair bonding behaviors, including latency to pair, and pairing stability 95,96 . In male cichlids, Amatitlania nigrofasciata, a IT/AVT antagonist cocktail inhibits affiliation with prospective partner and aggression towards nonpartners 43 . However, nonapeptides do not appear to be involved in male A. nigrofasciata pair bond maintenance 43,44 . Pair bonding pine voles, M. pinetorum, exhibit higher NAcc OTR expression 97 and VP V1aR densities 98 than do non-pairing montane voles, M. montanus. In five species of zebra finches (f: Estrildidae), LS V1aR density predicts species-typical social group sizes 51 . Similarly, in seven species of butterflyfishes (f: Chaetodontidae) AVT-ir neuron fibre varicosity density within the lateral septum-like region (the ventral and lateral parts of the ventral telencephalon, Vv/Vl) predicts species-typical pairing from non-pairing sociality 42 . Finally, in zebra finches, T. guttata, during pair bond formation and in established pairs, DA neurons expressing immediate early gene Fos (a marker of neuron activity) in the VTA is heightened 99,100 and pair bonded birds exhibit higher levels of DA in the ventral medial striatum (the super-structure of the NAcc) than unpaired birds 100 . However, since comprehensive examination of both the functional involvement of nonapeptide, dopaminergic, and opioid systems, and their respective targeting sites is thus far limited to M. ochrogaster, it is currently not possible to confidently discern the extent to which pair bonding regulatory neural networks may have converged across vertebrates.
The aim of this study was to determine whether IT, V1a, DA, and MO receptors are involved in pair bonding in a teleost, as well as establish the specific brain region(s) (anatomical substrate(s)) upon which each receptor type operates. Chaetodon lunutaus, was selected as a study species because it exhibits both pair bonding and solitary phenotypes among sympatric individuals of otherwise similar ecologies, offering an opportunity for comparative research. Receptor gene expression within eight distinct brain regions ( Table 1) was contrasted between pair bonded and solitary (control) individuals in order to reveal mechanistic correlates of pair bonding. Importantly, these regions include the putative ancestral homologs of those involved in M. ochrogaster pair bonding ( Table 1). Figure 1. The neurochemical and -anatomical substrates of M. ochrogaster pair bonding pre-date the split between ray-and lobe-finned fishes, and have since remained highly conserved. Hence, the convergence of vertebrate pair bonding may have relied on repeatedly co-opting these homologies that were already established in a common ancestor ~ 450 MYA. The distribution of vertebrate pair bonding is indicated by its presence (+) or absence (-) across major groups. The evolutionary origins and history of neurochemical systems (colored circles) and brain regions comprising the SDM (black circle) is shown. Abbreviations: AVT = arginine vasotocin; DA = dopamine; OP = opioid; IT = isotocin; SDM =social decision making network; MT = mesotocin; AVP = arginine vasopressin; OT = oxytocin; SDM = social decision making network. References for neural components: 1 , 1978). **In most mammals, the caudate putamen is a sub-structure of the striatum (O'Connell and Hofmann, 2011). ***While a tentative consensus for putative partial homologies between mammals and teleost brain regions has emerged, homologies should still be considered debatable Goodson and Kingsbury, 2013). ****Brain regions involved in M. ochrogaster pair bonding that were not examined in the current study include the PFC and the VP (because their ancestral homologs are unknown (O'Connell and Hofmann, 2011)), and the anterior hypothalamus (teleost homologue = vTn) (O'Connell and Hofmann, 2011).

Animal collection and sexing
This study was conducted at Lizard Island, located in northern section of the Great Barrier Reef (GBR), Australia (14°40′08″S; 145°27′34″E). To compare receptor gene expression within brain regions between pair bonded and solitary individuals, we first collected solitary and paired individuals of C. lunulatus from fringing reefs around Lizard Island. The social system of individuals was recorded following 5-min observations prior to collecting fishes by spearing through the dorsal musculature. Individual fishes were immediately placed in an icewater slurry for 5 min after which the brain was dissected (within 10 minutes of capture), embedded in optimal cutting compound (OCT), and frozen in liquid nitrogen for transportation to the laboratory where they were then transferred to -80 ͦ C freezer until sectioning. In order to sex individuals, gonads were removed and fixed in formaldehyde-acetic acid-calcium chloride (FACC) for at least one week. Thereafter, gonads were dehydrated in a graded alcohol series, cleared in xylene, embedded in paraplast, sectioned transversely (7 µm thick), and stained with hematoxylin and eosin. Sections were examined under a compound microscope (400 X magnification) for the presence of sperm (male) or oocytes (female) 101 .

Brain region extraction and measuring gene expression
Frozen brains were transversely sectioned on a cryostat at 110µm, thaw mounted onto Superfrost Plus slides (Fisher Scientific) and stored at -80˚C prior to brain region extraction (approx. one week). Brain regions, identified using a butterflyfish brain atlas 102,103 , were manually extracted at -30˚C using a hand-held micro-punching device (50mm diameter; Stoelting, model # 57401) (Fig. 2), incubated in RNAlater® at 4˚C over night, and then stored at -20˚C for up to one week. Brain region punching regime was standardized across individuals (see Supplimentary Table 1 for regime). Tissue punches were immediately transferred into a lysis buffer and homogenized by passing through a 21-gauge needle 15 times. Following, RNA was extracted using an E.Z.N.A® HP Total RNA kit (Omega, model # R6812-02) according to the manufacturer's instructions, and stored at -80˚C prior to cDNA synthesis. RNA was reverse transcribed into cDNA using Superscript III reverse transcriptase (Life Technologies) and genespecific primers (see Supplimentary Table 2 for primer sequences). Residual primers and salts from reverse transcription were removed using an E.Z.N.A ® Tissue DNA purification kit (Omega, product # D3396-02), and cDNA was stored at -20˚C for up to 10 days prior to qPCR.
The whole brain transcriptome of C. lunulatus was sequenced in order to use as a reference for designing species specific cloning primers for each gene of interest. One C. lunulatus brain was taken out of RNAlater, rinsed in 1X phosphate buffered saline (PBS) and placed immediately in Trizol (Life Technologies, Grand Island, NY) where RNA was extracted according to manufacturer instructions. Poly-adenylated RNA was isolated from each sample using the NEXTflex PolyA Bead kit (Bioo Scientific, Austin, TX, USA). Lack of contaminating ribosomal RNA was confirmed using the Agilent 2100 Bioanalyzer. A strand specific library was prepared using the dUTP NEXTflex RNAseq kit (Bioo Scientific), which includes a magnetic bead-based size selection of roughly 350 bp. The library was pooled in equimolar amounts with sample from an unrelated study after library quantification using both quantitative PCR with the KAPA Library Quantification Kit (KAPA Biosystems, Wilmington, MA, USA) and the fluorometric Qubit dsDNA high sensitivity assay kit (Life Technologies), both according to manufacturer instructions. Libraries were sequenced on an Illumina HiSeq 2000 to obtain paired-end 100bp reads. We first corrected errors in the Illumina reads using Rcorrector (parameters: run_rcorrector.pl -k 31) and then applied quality and adaptor trimming using Trim Galore! (http://www.bioinformatics.babraham.ac.uk/projects/trim_galore/; parameters: trim_galore --paired --phred33 --length 36 -q 5 --stringency 5 --illumina -e 0.1). After filtering and trimming, a total of 64,795,096 paired reads remained for de novo assembly. We created a C. lunulatus de novo transcriptome assembly using Trinity (parameters: --seqType fq --SS_lib_type RF). The raw Trinity assembly produced 376,338 contigs (N50: 1148 bp).
Using the C. lunulatus transcriptome as a reference, species-specific primers were designed to clone target gene sequences. Cloned target sequences were examined to determine whether they contained an exon-exon boundary using Danio rerio and Stegastes partitus complete genomes as a reference. If target gene sequences did not flank an exon boundary, then a second set of primers were designed to extend the obtained sequence towards the exon(s). Exon-containing sequences of ITR, V1aR, D1R, D2R and MOR genes were then used to design qPCR primers that flanked exon boundaries (18S ribosomal notwithstanding, because it does not contain an exon boundary) (see Supplimenary Table 2 for primer sequences). Prior to qPCR, primer sets and instrument cycling parameters were empirically optimized on standard curves using several metrics of quality control (i.e., assay amplification R 2 value of at least 0.95, assay slope of approximately -3.3, assay melting curve that only produced a single amplicon peak, no amplicon signal in the no template control (NTC) or nor reverse-transcription control (NRTC)). Quantitative PCR was then performed on each sample using a reaction mixture and qPCR cycling instrument (CFX380) that was recommended by the enzyme manufacturer (see Supplimenary Table 3 for parameters). Samples were run in technical triplicate on 384 well qPCR plates with standard curves in order to determine assay efficiency from the slope. Since assay efficiency was not the same across individual assays, averaged gene expression (C t ) values were standardized to assay efficiency prior to normalizing to 18S ribosomal RNA following methods of 104 . Not all focal regions of each brain were measured for gene expression due to insufficient tissue available.

Statistical analysis
For each gene within each brain region, a 2-way ANOVA with social system and sex as fixed factors was used to compare gene expression among treatment levels. Prior to analysis, data was natural log transformed +1 to improve normality of residual variance. Statistical analysis was conducted using SPSS software.

Nonapeptide receptors (ITR and V1aR)
Both ITR and V1aR expression within the Vv/vl differed interactively between sex and social system (Figure 3A, B; Supplimentary Table 4). In females, ITR and V1aR Vv/Vl gene expression was higher in pair bonding than solitary individuals (F 1,13 = 9.06, p< 0.01; F 1,13 = 9.18, p = 0.01, respectively); however, in males, there was no difference in either ITR or V1aR Vv/Vl gene expression between social systems (F 1,13 = .002, p = 1; F 1,13 = .036, p = 1, respectively). ITR and V1aR Vv/vl gene expression differed significantly between sexes (p < 0.05), with females having higher nonapeptide gene expression than males; and they also differed between social system (p < 0.05), with pair bonded individuals displaying higher nonapeptide gene expression than singletons. Gene expression of both ITR and V1aR did not differ between sex or social system interactively or independently in any other brain region, namely in the Dl, Dm, POA, TPp, Vc, Vc, or Vs (Supplimentary Table 4). V1aR gene expression was detected in all brain regions examined (ie., the Dl, Dm, POA, TPp, Vc, Vd, Vv/Vl, and Vs) and ITR gene expression was detected within all brain regions except for the Vd.

Dopamine receptors (D1R and D2R)
Gene expression of neither of the dopamine receptor class differed interactively between sex and social system in any brain region (Supplimentary Table  4). Furthermore, gene expression of dopamine receptor class did not vary between sexes (Supplimentary Table 4). However, in several brain regions, gene expression of both dopamine receptor classes differed between social systems (p <0.05 for each region), with both male and female pair bonded individuals expressing less than their solitary counterparts in these areas: D1R: Dm (trend), Dl, Vs, POA; D2R: Dm, Dl, Vs, POA, Vc, and TPp (Fig. 4A, B; Supplimentary Table 4). Dopamine receptor class expression did not vary between social systems in other brain regions: D1R: TPp, Vd, Vv/vl; D2R: Vd, Vv/vl (Supplimentary Table 4). D1R and D2R gene expression was found in all brain regions examined (ie., the Dl, Dm, POA, TPp, Vc, Vd, Vv/Vl, and Vs). . ITR (A) and V1aR (B) gene expression differences between sexes and social systems of C. lunulatus within the ventral and lateral parts of the ventral telencephalon (Vv/Vl). Boxes show the first and third quartiles, the black line in box represents median, and whiskers represent minimum and maximum value. Sample sizes are listed below each treatment group. Asterisks indicate statistically significant differences between treatment groups (P < 0.05) 2-way ANOVA and HSD Tukey Test.

Mu-opioid receptor (MOR)
Gene expression of MOR did not differ interactively between sex and social system in any brain region, nor did it differ independently between sexes in any brain region (Supplementary Table 4). However, in several brain regions, MOR gene expression differed between social systems (p < 0.05 for each region), with both male and female pair bonded individuals expressing less than their solitary counterparts in the Dl, Vs, POA, and TPp (Fig.  5A, B; Supplementary Table 4). MOR receptor expression within the Vc trended lower in pair bonded individuals than in solitary counterparts; however, this was to a statistically nonsignificant extent (p = 0.059). MOR receptor expression within the Dm, Vd, and Vv/Vl did not differ between social systems (Supplementary Table 4). MOR gene expression was detected in all brain regions examined (ie., the Dl, Dm, POA, TPp, Vc, Vd, Vv/Vl, and Vs).

Nonapeptide circuitries of pair bonding in Chaetodon lunulatus
Our comparative analyses reveled that in females, paired individuals displayed higher ITR and V1aR receptor expression within the Vv/Vl than solitary individuals, indicating that ITR and V1aR signaling within the Vv/Vl might be important for mediating pair bonding for this sex. However, in males, receptor gene expression did not differ within any brain region between social conditions. Few other studies have explored the involvement of nonapeptides in teleost pair bonding, and they have been mostly on males. In male cichlids, A. nigrofasciata, a general nonapeptide receptor antagonist inhibits affiliation with a prospective partner and aggression towards non-partners 43 , indicating the involvement both systems in pair bond formation. However, nonapeptide signaling does not appear to be involved in pair bond maintenance in males of this species 43,44 . In established pairs of Neolamprologus pulcher (but not of Telmatochromis temporalis) cichlids, whole brain gene expression of IT is positively correlated with partner affiliation 105 . Additionally, N. pulcher displays higher whole brain gene expression IT than the less affiliative Telmatochromis temporalis 105 . Similar to our study, Dewan et al. (2011) found in males of seven species of chaetodontids, Vv/Vl AVT-ir neuron fibre varicosity density predicts species-typical pairing from non-pairing sociality. Taken together, these studies indicate that while nonapeptides play a recurring role in promoting teleost pair bonding, this is species-, gender-, and context-specific.
What might be the precise functional role of Vv/Vl IT-ITR and AVT-V1aR signaling in promoting female C. lunulatus pair bonding? In teleosts, both AVT and IT mediate a wide range of behavioral domains that that lack a universal valence and appear to be context specific 106 . AVT regulation of social behavior has been studied extensively, albeit almost exclusively in males, where it has functionally been shown to promote spawning [107][108][109] , mate guarding or other forms of conspecific aggression/avoidance 43,108,[110][111][112][113][114][115] , social communication 116 , social preference 117 , and approach behavior 118 . Few studies have functionally examined IT involvement in teleost social behavior and again, those which have are exclusive to males. Similar to AVT, these studies demonstrate an inconsistent effect of IT on approach and avoidance behaviors, including pair bond formation and maintenance 43,44 , social preference 117 , and social affiliation 119 . However, since neurochemical effects are often specific to the site(s) of action 120,121 and sites are not confirmed in these studies; it is difficult to know which, if any, of these roles generate insight into the current findings. The rich literature on rodents, however, shows that in pair bonding species, AVP within the lateral septum (LS, the mammalian homolog of Vv/Vl) is involved in both partner affiliation 24 and territoriality 122 , perhaps reflecting its broader role in social recognition/memory [123][124][125][126][127][128] . As with AVP, septal (including lateral septal) OT is essential for social recognition in rodents [128][129][130][131] . In pair bonding butterflyfish, both olfactory and visual cues are used for conspecific recognition 132 and are necessary to modulate relationships with partners, territorial intruders 133 and competitors for mates 134 . In teleosts, the ventral telencephalon is the major target of olfactory projections [135][136][137] , but not of optic neurons, nor is it innervated with the optic tectum [138][139][140] , the major brain region in which visual information is integrated and processed in vertebrates 141 . Taken together, we speculate that in C. lunulatus females, V1aR and ITR activation within the Vv/Vl might enhance conspecific recognition via olfactory perception 42 . This is certainly an intriguing area of further inquiry.

Convergent evolution with birds and mammals
The teleost, mammalian, and avian lineages share striking similarities in nonapeptidemediated pair bonding circuitry. We have shown here that, as in other teleosts 42,43 , AVT plays an important role in C. lunulatus pair bonding, and that its effects are likely exerted within the Vv/Vl through V1aR activation, mirroring the role of AVP in birds 142 and of AVP-V1aR binding within the LS (the mammalian homolog of Vv/Vl) in M. ochrogaster voles 24 . Similarly, fMRI studies show that in humans, activation of the septum, which is rich in AVP binding sites 143 is associated with "obsessive love" 144 . We have further discovered that, as in other teleosts, IT is important for C. lunulatus pair bonding, paralleling the functional involvement of OT in pair bonding M. ochrogaster rodents 21,23,145 and in non-human primates (i.e., marmosets, Callithrix penicillata 93 ; tamarins, Saquinus oedipus 94 ).
As with AVP, OT activity is also implicated in human pair bonding: intranasal OT in men within romantic partnerships increases preferred interpersonal distance from non-partner females 146,147 , and plasma OT levels predict future success rates in romantic relationships 148 . However, unlike the AVP/AVT system, the site(s) of OT action in humans and M. ochrogaster rodents (ie. the prefrontal cortex (PFC), and nucleus accumbens (NAcc) (the mammalian homologue of the Vd) 21,149,150 are different than that of IT in teleosts (i.e., the lateral septumlike area). There are several potential explanations for this. First, since the evolutionary antecedent of the mammalian PFC is unclear 151,152 , it couldn't be examined here. Secondly, and somewhat surprisingly, this study found no ITR expression in the NAcc/Vd at all, suggesting that this region is not important for pair bonding and social behavior in general in C. lunulatus. Alternatively, this could be an artifact of lack of tissue available for sampling. Given that ITR within the NAcc/Vd is expressed in teleosts 153 and that the NAcc/Vd is considered a key node in the vertebrate SDM network, we suspect that the latter alternative is more plausible. Technical limitations notwithstanding, we might still expect anatomical targets of ITR/OTR-mediated pair bonding to be distinct between mammals and teleosts, due to differences in pre-existing neural circuitries that would have been available for co-option during their independent evolution. In mammals, a pre-established OT-mediated maternal bonding circuitry, in which the NAcc is a critical site of action 154,155 is thought to have been recruited during the evolution of pair bonding [156][157][158][159] . Since parental care did not precede the evolution of pair bonding in butterflyfishes 134 , this pre-existing circuitry would have been unavailable for co-option in these organisms.

Dopaminergic circuitries of pair bonding in Chaetodon lunulatus
An essential component of pair bonding is the reinforcement of partner affiliation 11 , which relies on individuals to perceive their partners as rewarding (i.e., approach eliciting) through heightened "salience" 32,160,161 . In pair bonding prairie voles, conspecific affiliation is not naturally rewarding, so does not facilitate pair bonding independently 11 . However, affiliation coupled with natural reward (specifically mating), is reinforcing and thus promotes pair bonding 162 . Therefore, mammalian pair bond formation is viewed to depend on conditioned reward learning, whereby individuals learn to associate their partner (conditioned stimulus) with mating (natural reward/unconditioned stimulus) 31,163-165 . The associative reward learning involved in this conditioned partner preference (CPP) is dependent upon dopamine acting upon nodes of the mesolimbic reward system--a neural network where the salience of environmental stimuli is evaluated 8,21,161 . Similar to pair bonding prairie voles, in situ partner removal experiments on C. lunulatus show that widowed males and females initially act antagonistically when approached by opposite sexed conspecifics within their territory. However, persistent "stalking" by the intruder towards the widowed individual while foraging accompanies the development of a new pair bond (Nowicki et al., manuscript submitted). Hence, C. lunulatus pair bonding might also rely on the learned association between partner (conditioned stimulus) and food (natural reward/unconditioned stimulus), and this associative learning might also be underpinned by dopamine acting upon reward circuitry. In support of this idea, in teleosts, both DA-D1R and -D2R binding are critical for psychostimulant/food reward learning 76,77,[166][167][168][169] and a network structured very similarly to the amniote mesolimbic reward system has been identified 88,89 . Importantly, almost all of the brain regions that were associated with DA-mediated pair bonding in this study are nodes of this putative teleost mesolimbic reward system (POA notwithstanding).
Our comparative results revealed that in both male and female C. lunulatus, D2R gene expression differed in the TPp and Vc between pair bonded and solitary individuals. This suggests that DA-D2R signaling within these regions may be important for pair bonding in both sexes of this species. The mammalian homologs of these brain regions, namely the ventral tegmental area (VTA) and striatum (STR), comprise the central ascending dopaminergic innervation pathway in the mesolimbic reward system 89,170 . In mammals, this pathway appears to have been co-opted during the evolution of pair bonding in order to mediate partner reward learning 11,30 . In teleosts, the TPp seems to have the densest cluster of DA-synthesizing cell bodies in the brain 170 ) and is considered the dopaminergic system ascending to the striatum 80 . Furthermore, DA-synthesizing neurons within the TPp are necessary for conditioned learning of place preference 77,171 . Given the aforementioned homologies and functional similarities, we tentatively hypothesize that DA ascending from the TPp and binding to D2Rs within the striatal Vc might function to mediate partner--consumatory reward learning in a similar manner to the VTA-NAcc complex in mammals. Yet the hypothesis that the TPp is a major source of DA in C. lunulatus pair bonding does not explain why it appears to be a potential target of DA action in this species? Perhaps the TPp is both a source and a site of action in DA-mediation of pair bonding in C. lunulatus. This possibility might also apply for mammalian counterparts, because the VTA-mPFC complex within the mesocorticolimbic pathway is reciprocally innervated [172][173][174] and DA-synthesizing neurons within the VTA display a high density of dendrite D2 receptors 64 .
Our comparative results revealed several other potential sites of DA action that are shared by D1R and D2R targeting. This is consistent with the ideas that D1R modulates D2R mediated events 175 and that D1-and D2R subtypes function complementarily to mediate pair bonding behavior 32 . Three of these implicated brain regions, the Dm, Dl, and Vs, belong to the putative teleost mesolimbic reward system, and mediate emotional learning/memory 81 , relational/spatial/temporal memory 81 and aggression/spawning 88 , respectively. The final brain region implicated, the POA, is a node of the conserved social decision making neural network, where it mediates several social domains across vertebrates, including sexual activity and male aggression 88,[176][177][178][179][180] . In voles, in particular, the mPOA appears critical for several pair bonding behaviors, including pair bond formation 181 , mating 180 , mate guarding, and territorial defense 177,178,182 . mPOA-mediated pair bond formation and mating in particular are believed to be attributed to dopamine 182 .
Hence, we propose that in C. lunulatus, D1 and D2 receptors might act synergistically within the Dm, Dl, Vs, and POA to mediate emotional, spatial/temporal, and sexual/mateguarding mnemonic events involved in partner reward learning. Of final note, dopamine receptor expression within these brain regions was relatively lower in pair bonded fish than in solitary counterparts. This is somewhat contradictory, because since we hypothesize that DA binding promotes partner reward learning, signaling is expected to be relatively higher in paired fish. We offer two potential explanations for this. First, reduced DAR gene expression might reflect reduced receptor expression, which might act as a compensatory mechanism for heightened DA release 183 . Secondly, while gene expression often increases with the activity or abundance of protein products, it has also been shown to exhibit an inverse relationship 184 , as has been previously shown in pair bonding M. ochrogaster 185 . Hence, it is possible that relatively lower DAR gene expression (and MOR gene expression, see below) reflects relatively higher levels of receptor abundance or activation in association with C. lunulatus pair bonding.

Convergent evolution with birds and mammals
The teleost, bird, and mammalian lineages share some prominent similarities in dopamine-mediated pair bonding circuitry. Our comparative results suggest that dopamine neurotransmission within the mesolimbic reward network is important for pair bonding in C. lunulatus, as appears to be the case in the zebra finch T. guttata 99,100,186,187 and in M. ochrogaster rodents [30][31][32] . A notable brain region of this network that appears to be targeted by DA in all three taxa is the striatal Vc/ striatal NAcc 31,32,100,186 . In addition, DA appears to act within the TPp (mammalian and avian VTA) and the POA in both C. lunulatus and T. guttata 99,186 , but whether it targets these regions in mammals remains untested. Finally, our study further implicated that DA-D1R-and -D2R signaling within the Dm, Dl, and Vs might also regulate C. lunulatus pair bonding, but their involvement within homologous regions (i.e., the blAMY, HIP, and meAMY/BNST, respectively) remain untested in other taxa. Interestingly, however, a growing body of research implicates that DA targets similar regions of the mesolimbic reward system to regulate partner affiliation in humans 188,189 . For example, functional magnetic resonance imaging (fMRI) shows that striatal regions, as well as the VTA, AMY, and the HIP, which are rich in dopamine activity 64 , are activated differently when participants view images of those with whom they're in an intense romantic or long-term, deeply-loving relationship; than when viewing pictures of other familiar individuals [190][191][192] .

Mu-opioid receptor circuitry of pair bonding in Chaetodon lunulatus
Mu-opioid receptor gene expression varied in relation to pairing sociality within the POA and several nodes of the mesolimbic reward system: the striatal Vc, Dl, Vs and TPp (however, in the striatal Vc this was to a statistically non-significant extent (p = 0.059). In teleosts, the POA mediates social and feeding behavior 88 . It is well established in mammals, that MOR plays an essential role in mediating the reinforcing effects of natural rewards (e.g., food, water, sex, social affiliation) and of psychostimulant rewards by eliciting motivational and pleasurable hedonic responses to these stimuli 68,[193][194][195][196][197][198][199][200][201][202][203][204][205] . Preliminary investigations suggest that opioid and/or MOR action within mesolimbic reward system also mediates reward processing in teleosts 76,77 .
In prairie voles in particular, the mu-opioid system plays a critical role in facilitating pair bond formation [35][36][37] . Its effects are exerted within striatal regions of the brain, including the dorsal striatum, dorsomedial NAcc shell, and caudate putamen [35][36][37] ; where dorsal striatum MORs are believed to facilitate pair bond formation by promoting mating during CPP, and dorsomedial MORs are believed to do so by modulating the positive hedonics of mating during CPP 36 . (See section on dopamine for description of CPP cognitive process.) In teleosts, food is a natural reward whose reinforcing properties are modulated by the opioid system 76 . In C. lunulatus, exclusive pair-wise feeding strongly coincides with pair bond formation and maintenance (Nowicki et al., manuscript submitted). Hence, we hypothesize that OP-MOR binding within the POA and nodes of the mesolimbic reward system (i.e., the Vc, Dl, Vs, and TPp) promotes pair bonding in C. lunulatus by modulating the positive hedonics of natural consumatory reward during CPP learning. In further support of this idea, our comparative results revealed that the opioid and the dopaminergic systems appear to target several of the same nodes of the mesolimbic reward system (i.e., the Vc, Dl, Vs, and TPp), indicating that they might converge on these regions in order to underpin the learned association between consumatory reward affect and one's partner, respectfully, during the CPP process. Experimental research is needed to empirically test this hypothesis.

Convergent evolution with mammals
To date, potential involvement the opioid-mu-opioid system in pair bonding has only been examined in two species, C. lunulatus butterflyfishes (current study), M. ochrogaster voles 33,[35][36][37]206 . In both organisms, it appears to play an important role, and effects seem to be exerted by acting upon nodes of the mesolimbic reward network. Specifically, we found comparative evidence that the striatal Vc -MORs are important in C. lunulatus, mirroring the role of striatal region (the NAcc and dorsal striatum) MORs in M. ochrogaster 36,37 . While several other nodes of the mesolimbic reward system (i.e., the Dl, Vs, TPp) and the POA were implicated in C. lunulatus pair bonding, the involvement of their homologs (i.e., the meAMY, LS, VTA and POA, respectively) in other taxa is yet to be explored functionally, so cannot be compared here. Interestingly, however, emerging evidence suggests that OP-MOR activity is important for pair bonding in humans as well, where here too its function is believed to be eliciting motivation and positive hedonics in response to romantic affiliation [207][208][209] . Moreover, the implicated brain regions involved share some similarities with those implicated in C. lunulatus and established in M. ochrogaster. Specifically, fMRI studies suggest that in humans, motivational aspects of partner preference formation are regulated by the dorsal striatum 206 , which is rich in MORs 210 . Whereas, the positive hedonics of "romantic love" are associated with the AMY, septal fornix, and VTA 190,192 , all of which are also rich in MORs 211,212 . Working model for the neural network of pair bonding in fishes By synthesizing our current findings with available information on teleost neurochemical synthesis and projection pathways, and functional insight from pair bonding M. ochrogaster counterparts, we can now begin to assemble a working model for how isotocin, arginine vasotocin, dopamine, and opioid systems might interplay to comprise a broader neural network of C. lunulatus pair bonding (as illustrated in Fig. 6). It is important to emphasize from the very onset that the only component of this working neural network model that derived from our findings is the involvement of IT-ITR, AVT-V1aR DA-DR, and OP-MOR signaling within brain regions (olfactory bulb notwithstanding), and that the remainder of this model is purely speculation.
We tentatively hypothesize that pair bonding in C. lunulatus relies on conditioned partner preference (CPP), as in the mammalian model, M. ochrogaster 11 . Several lines of behavioral evidence support this hypothesis. Field observations reveal that solitary C. lunulatus do not find prospective partners naturally rewarding (i.e., they respond to prospective partners by antagonism rather than approach). Only after continued and exclusive cohabitation involving pair-wise foraging, is a pair bond developed (Nowicki et al., manuscript submitted). After development, pair bonds are enduring, and are characterized by selective affiliation and feeding with partner, and selective antagonism towards nonpartners (Nowicki et al., manuscript submitted). We propose that during CPP in C. lunulatus, individuals form a learned association between natural reward (food, unconditioned stimulus) and their new partner (conditioned stimulus) resulting in the new partner to take on rewarding properties, and thus reinforce selective approach behavior 8 . In our working model for this process, feeding activates the TPp (VTA), concurrently triggering OP-MOR and DA-D2R activity within the mesocorticolimbic reward system, which converges in the Vc (striatum), Dl (HIP), Vs (meAMY, BNST), TPp (VTA), and POA in particular, thereby modulating consumatory reward affect and reward learning/salience of partner-associated cues. In this pathway, the major source of DA projection to at least the Vc (striatum) is most likely the TPp (VTA) 44,64 , while OP is most likely to originate from the hypothalamic nucleus lateralis tuberis 213 . Meanwhile, olfactory cues from the partner are transmitted from the olfactory bulb (OB), ultimately reaching the Vv/Vl (LS), where IT and AVT nonapeptide activity converges to promote olfactory learning in females. The source of nonapeptide release originates from cell bodies within the POA [214][215][216] . After pair bond formation, concordant D1R and D2R activity within the mesocorticolimbic reward system and POA modulates pair bond maintenance by mediating aggression towards non-partner conspecifics 206 . Also involved in this neural network would be higher-order motor circuits that underpin the behavioral outcome of approach and affiliation towards partner, and aversion towards non-partners (not studied nor illustrated here) 217 . This proposed working model is speculative, and requires empirical testing (see below). Figure 6. Sagittal view of brain illustrating a working neural network model for pair bonding in Chaetodon lunulatus. Colors within brain regions represent putative sites of action for each system based on comparative data on receptor gene expression provided here. Symbols and arrows indicate the sex(es) to which receptor phenotypes apply, and direction of receptor gene expression (up-or down-regulated) within brain regions, respectively. Colored lines represent putative neurochemical projections from predominant sites of synthesis (based on literature) to putative target sites (based on current findings). Illustration made by J.P.N., adapted from Dewan and Tricas, 2014, with permission. This is one of the first studies to explore the neurobiology of pair bonding in an early vertebrate (i.e., a reptile or an anamniote). While few other studies have researched the involvement of nonapeptides (i.e., teleosts: AVT: 42,43,105,218 ; IT: 43,44 ); that we are aware of, this is the first to research the involvement of the dopamine and opioid systems, and examine gene expression in specific brain regions. Although we provide support for the involvement of these neurochemicals, their conjugate receptors, and brain regions; our data are only correlative. Furthermore, due to the paucity in neural research on fish pair bonding, we have relied heavily on drawing upon the rich body of literature on the mammalian model, Microtus ochrogaster, in order to speculate the cognitive and behavioral functions that these putative neural circuits might subserve. Therefore, it is important to consider that our proposed neural network model of teleost pair bonding is far from conclusive and is certainly incomplete. Nonetheless, we believe it provides a useful foundation from which specific hypotheses related to teleost pair bonding can be tested in the future. A priority should now be to experimentally validate whether the neuroanatomical correlates of pair bonding found here are functionally relevant, and if so, then whether these functions are analogous to those of mammalian counterparts. Furthermore, we advocate exploring the potential involvement of other promising brain regions that are critical to vertebrate social behavior, including the periaqueductal gray/central gray (PAG/CG), ventral tuberal nucleus (vTn) (homologous to the mammalian anterior hypothalamus, (AH)) and anterior tuberal nucleus (aTn) (homologous to the mammalian ventromedial hypothalamus, (VMH)) 88,91,92 . Moreover, similar preliminary investigations into the involvement of other likely candidate systems that modulate reward and positive reinforcement behavior (e.g., serotonin 10 and orexin 219 ), negative reinforcement behavior (i.e., corticotrophin releasing factor 220,221 ), and motor output (e.g., GABAergic and glutamatergic) 217 ) are encouraged. Finally, in order to better understand the extent to which neurobiological systems of pair bonding have converged across vertebrate evolution, complementary research needs to be done on multiple species within and across all major taxonomic groups.

Conclusions
In addition to representing an integral part of the human experience [3][4][5]8 , pair bonding has independently evolved in every major vertebrate lineage. It is already clear that this has not occurred through a compete and universal convergence of a single regulatory neural network. However, our study contributes to an emerging pattern that in at least selective cases, even those involving phylogenetically distant taxa with distinct evolutionary histories, this might occur through at least a partially converged neural network. M. ochrogaster rodents and C. lunulatus teleosts, despite being separated by ~450 million years of independent evolution 47 , and despite having opposing parental evolutionary histories (parental vs. non-parental, respectively), appear to share some striking similarities in the neural substrates that underpin pair bonding. We have discovered evidence for the involvement of isotocin, arginine vasotocin, dopamine, and opioid systems in C. lunulatus pair bonding, corresponding to their (or their homologs) involvement in M. ochrogaster counterparts. Moreover, we have described that in association with pair bonding sociality, nonapeptide receptor expression varies in the lateral septum-like region, while dopamine and opioid receptor expression varies within other regions of the mesolimbic reward network, including the striatum; mirroring sites of action in M. ochrogaster. Therefore, we tentatively suggest that the neurobiology of pair bonding between these taxa has at least partially converged through the repeated co-option of evolutionarily deep molecular and anatomical homologies that were already established in ancestral osteichthyes ~ 450 MYA. In order to determine the extent to which this has occurred across vertebrates, complementary studies across a wider range of lineages (most urgently amphibians and reptiles) are now needed.