Transcriptional regulatory divergence underpinning species-specific learned vocalization in songbirds.

Learning of most motor skills is constrained in a species-specific manner. However, the proximate mechanisms underlying species-specific learned behaviors remain poorly understood. Songbirds acquire species-specific songs through learning, which is hypothesized to depend on species-specific patterns of gene expression in functionally specialized brain regions for vocal learning and production, called song nuclei. Here, we leveraged two closely related songbird species, zebra finch, owl finch, and their interspecific first-generation (F1) hybrids, to relate transcriptional regulatory divergence between species with the production of species-specific songs. We quantified genome-wide gene expression in both species and compared this with allele-specific expression in F1 hybrids to identify genes whose expression in song nuclei is regulated by species divergence in either cis- or trans-regulation. We found that divergence in transcriptional regulation altered the expression of approximately 10% of total transcribed genes and was linked to differential gene expression between the two species. Furthermore, trans-regulatory changes were more prevalent than cis-regulatory and were associated with synaptic formation and transmission in song nucleus RA, the avian analog of the mammalian laryngeal motor cortex. We identified brain-derived neurotrophic factor (BDNF) as an upstream mediator of trans-regulated genes in RA, with a significant correlation between individual variation in BDNF expression level and species-specific song phenotypes in F1 hybrids. This was supported by the fact that the pharmacological overactivation of BDNF receptors altered the expression of its trans-regulated genes in the RA, thus disrupting the learned song structures of adult zebra finch songs at the acoustic and sequence levels. These results demonstrate functional neurogenetic associations between divergence in region-specific transcriptional regulation and species-specific learned behaviors.

Introduction

Species-specific behavior plays a role in a variety of inter- and intraspecific interactions, including reproduction and habitat use, in which species differences are thought to be an important factor in species co-occurrence [1-3]. Such species-specific behaviors can arise via species differences in the structure and development of the neural circuits underlying behavior [4-6]. Differences between closely related species are thought to be driven by differential expression and functional changes of orthologous genes in conserved neural circuits, which are often in turn driven by transcriptional regulatory divergence [7-10]. Transcriptional regulatory divergence between species can arise due to species divergence in cis-regulatory elements that affect the transcriptional rate and stability, and/or in trans-regulatory factors that access cis-regulatory elements [11-16] (). However, it remains largely unknown how transcriptional regulatory divergence contributes to the generation of species-specific behavior, especially in the case of learned behavior.

Songs produced by oscine birds are complex vocal signals acquired through vocal learning [17,18]. Songs are species-specific, and these species differences play an important role in mating interactions and territory defenses within and between species [1,19,20]. In the songbird brain, a conserved neural circuit specialized for vocal learning, called the song system, contributes to song learning and production [18,21,22]. Birdsong is composed of two main traits associated with species specificity: the acoustic elements (syllables) and the temporal pattern (sequence) of song. The production of syllable acoustics and sequence is mainly regulated by the robust nucleus of the arcopallium (RA) and the song nuclei HVC (proper name), respectively, in the vocal motor circuit of the song system () [22-24]. The importance of these song nuclei in determining species-specific song traits suggests an underlying causative role of species differences in the structure and activity of these regions. Consistent with this, a variety of genes, including transcription factors and neuromodulator receptors, are differentially expressed in these song nuclei between species, even in a laboratory-controlled environment [25-27]. However, a key gap in our knowledge is how species-specific patterns of gene expression in these regions arise via regulatory differences between species.

In this study, we used two closely related songbird species, zebra finch (ZF; Taeniopygia guttata), owl finch (OF; T. bichenovii), and their interspecific first-generation (F1) hybrids, to elucidate how transcriptional regulatory divergence is associated with species-specific song (). These two species diverged about 6.5 million years ago and share overlapping habitats in the north and west of Australia [28,29]. In addition, they produce songs with characteristic species-specific syllable acoustics and sequence. By comparing the gene expression ratio between the two species and the allele-specific expression (ASE) ratio in the F1 hybrids (), we assessed the total number of genes whose expression differs by divergence in cis- versus trans-transcriptional regulation between the two species. On the basis of Gene Ontology (GO) enrichment and the upstream regulatory analyses of transcriptional regulation–altered genes, we identified the candidate key upstream modulators of these differentially regulated genes and examined the functional effects of altered transcriptional regulation in the song nuclei.

Results

Species difference in song phenotypes between ZF and OF

First, we compared the song features of ZF and OF reared with conspecific song tutoring in our breeding colony to confirm whether a laboratory-controlled environment could maintain species-specific song features. We compared the songs of the two species regarding syllable acoustics and sequential features (12 parameters) at the adult stage () and identified significant differences in six acoustic syllable parameters (i.e., syllable duration, inter-syllable gap duration, entropy variance, amplitude modulation [AM] variance, mean frequency modulation [FM], and FM variance) and in syllable sequence features (motif and repetition transition rates) (n = 6 birds each, p < 0.01, one-way ANOVA) ( and ) [30, 31]. We found that the range but not the pattern of each acoustic feature’s distribution overlapped between ZFs and OFs (3,000 syllables from n = 6 birds each and 500 syllables/bird) (), thus suggesting that the species differences in the syllable acoustics were not caused by physical species-specific constraints in the peripheral vocal organs.

Species difference in song structures between ZF and OF.

(A) Typical examples of songs from ZFs and OFs that were reared with conspecific song tutoring and cross-species song tutoring. Orange solid and blue dotted lines represent the motif and repetitive structure of syllables, respectively. (B) Species differences in the syllable sequence of ZF and OF songs. (Left) Syllable similarity matrices for songs produced by ZFs and OFs that were reared with conspecific song tutoring and cross-species song tutoring. (Right) Motif and repetition indices of ZF and OF songs (n = 6 each from conspecific song tutored ZF and OF, n = 4 and 3 from cross-species song tutored ZF and OF, respectively; one-way ANOVA, *p < 0.05, **p < 0.01). Each dot corresponds to an individual bird. (C) Species differences in syllable acoustics (syllable duration, inter-syllable gap duration, entropy variance, AM variance, mean FM, and FM variance) of ZF and OF songs (“Con”: n = 6 each from conspecific song tutored ZF and OF; “Cross”: n = 4 and 3 from cross-species song–tutored ZF and OF, respectively; one-way ANOVA, *p < 0.05, **p < 0.01, ***p < 0.001). Each dot corresponds to an individual bird. Relevant data values are included in . AM, amplitude modulation; FM, frequency modulation; OF, owl finch; ZF, zebra finch.

We further performed cross-species song tutoring experiments to examine how genetic and environmental factors contribute to generate species-specific song features of ZF and OF (). Under the cross-species song tutoring condition, juveniles heard only the counter-species songs as tutor songs. By comparing the p-values of song feature differences between conspecific and cross-species song tutoring conditions, we found that song tutoring affected most of the song parameters, including syllable sequence and acoustics (i.e., inter-syllable gap duration, entropy variance, AM variance, mean FM, and FM variance) (). However, except for AM variance, all song parameters retained species specificity (ZF, n = 4, OF, n = 3; one-way ANOVA, p < 0.05). In line with this result, we performed principal component analysis (PCA) to investigate the song feature distribution of conspecific and cross-species song tutored birds by reducing the dimensionality of the syllable acoustics and sequential features. We observed that clusters were separable by species but not by song tutoring conditions (). As many studies in songbirds reported [32-35], these results also indicate that song learning of these two species is implemented based on species-specific genetic constraint.

Genome-wide transcriptional analysis between ZF, OF, and F1 hybrids

We then conducted a genome-wide transcriptional analysis to elucidate divergence of transcriptional regulation between ZF and OF in their song nuclei. For this purpose, using laser microdissected HVC and RA tissues from ZFs and OFs, we identified 11,501 and 11,487 genes in HVC and RA, respectively, as genes with detectable expression levels in either ZF or OF (reads per kilobase of transcript per million reads mapped [RPKM] ≥ 1). We then calculated the expression ratio between ZF and OF for each gene as “A/B” = RPKM(ZF average)/RPKM(OF average) (n = 4 birds each) (Figs ).

Species differences in gene expression between ZF and OF and ASE in F1 hybrids.

(A) Examples of species differences in gene expression between ZF and OF. (Left panels) Expression levels of GRIK1, RAB5A, and LMO7 in song nucleus RA of ZFs and OFs. Gray boxes represent the position of exons for each gene. Dark blue peaks below exons represent read density. (Right panels) Gene expression levels in ZF and OF and the average of the expression ratio between ZF and OF. Each dot represents the RPKM value for individual. Mean ± SEM (n = 4 birds each; one-way ANOVA, *p < 0.05, ***p < 0.001, n.s., not significant). RAB5A is an example with no expression difference between ZF and OF. (B) Examples of ASE in F1 hybrids. (Upper panels) Allelic expression ratios in F1 hybrids at species-specific SNPs (ss-SNPs) of RASGEF1B and HTR1B in song nucleus RA. Dark blue peaks below exons represent read density. White bars in the dark blue–colored peaks represent ss-SNP positions. Pie charts of each ss-SNP represent the percentage of transcribed read numbers from ZF (orange) and OF (blue) alleles. (Bottom panels) The percentage and ratio of parental species-allelic expression of RASGEF1B and HTR1B in OZ and ZO F1 hybrids. Each dot represents average allelic expression ratios of all ss-SNPs in one individual (n = 4 birds each, mean). Orange- and blue-colored bars represent the values from ZF and OF alleles, respectively. Mean ± SEM (n = 4 birds each). Relevant data values are included in . ASE, allele-specific expression; Chr, chromosome; F1, first-generation; GRIK1, Glutamate receptor, ionotropic, kainate type 1; HTR1B, 5-hydroxytryptamine receptor 1B; LMO7, LIM domain only protein 7; OF, owl finch; OZ, F1 hybrid offspring between OF♀ and ZF♂; RA, robust nucleus of the arcopallium; RAB5A, Ras-related protein Rab5A; RASGEF1B, Ras-GEF domain-containing family 1B; RPKM, reads per kilobase of transcript per million reads mapped; ZF, zebra finch; ZO, F1 hybrid offspring between ZF♀ and OF♂.

Based on a comparison of whole brain transcriptome between ZF and OF, a total of 2,409,063 SNPs were identified as species-specific SNPs (ss-SNPs) in their transcribed sequences. Using the ss-SNPs for the quantification of ASE ratios in the F1 hybrids, we set a cutoff to extract genes with ≥5 reads at each ss-SNP position and median RPKM ≥ 10 (n = 4 each from ZO and OZ hybrids). Totals of 5,827 and 6,328 genes passed the criteria in HVC and RA, respectively. The ASE ratio of each gene in individual F1 hybrids was calculated as “a/b” = Reads(ZF allele)/Reads(OF allele) (Figs ). To date, there is no evidence for paternal and maternal genomic imprinting in avian species [36]. In line with this, we identified no genes with a significant paternal or maternal bias in allelic expression in ZO and OZ hybrids. Furthermore, the two reciprocal F1 hybrids (ZO and OZ) have an extremely high correlation in their ASE ratios (Pearson correlation coefficient, r = 0.527, p < 2.2 × 10−16 in HVC; r = 0.550, p < 2.2 × 10−16 in RA) (). Therefore, we treated ZO and OZ hybrids equally when calculating ASE ratios.

Transcriptional regulatory divergence between ZF and OF

Transcriptional differences, cis- and/or trans-regulation, for each gene can be evaluated using the gene expression ratio between two species and the ASE ratio in the F1 hybrids [14,16,37-39]. ASE in the F1 hybrids reflects cis-dependent differences between the alleles of each parental species, because the two alleles of each gene are exposed to same trans-acting regulatory environment in cells. By comparing the gene expression ratio between parental species and the ASE ratio in F1 hybrids, we determined the following five categories of transcriptional regulatory divergences: (i) “cis-regulation” for genes with significant cis- but not trans-effects (with a/b ≠ 1 and A/B = a/b) as “cis-regulated genes,” (ii) “trans-regulation” for genes with significant trans- but not cis-effects (with a/b = 1 and A/B ≠ a/b) as “trans-regulated genes,” (iii) “both cis- and trans-regulation” for genes with both significant cis- and trans-effects (with a/b ≠ 1 and A/B ≠ a/b) as “both cis- and trans-regulated genes,” (iv) “conserved regulation” for genes with no significant cis- or trans-effects (with a/b = 1 and A/B = a/b) as “conserved genes,” and (v) ambiguous regulation (Figs ). For this categorization of transcriptional regulatory divergence, we applied a cross-replicate comparison of ASE ratios in the F1 hybrids, through which we could minimize incorrect estimation of cis- and artificial negative correlation in cis- versus trans-comparison (see Materials and methods) [40]. After this procedure, we observed that over 75% and 10% of the examined genes were expressed in both HVC and RA with either “conserved” or “ambiguous” regulation between ZF and OF, respectively (). In contrast, transcriptional regulatory divergence changed the expression of 158 (2.7% of the total 5,827 genes), 271 (4.7%), and 183 (3.1%) genes in HVC categorized as cis-, trans-, and both cis- and trans-regulated genes, respectively. Likewise, in RA, the expression of 246 (3.9% of the total 6,328 genes), 383 (6.1%), and 183 (2.9%) genes was altered by cis-, trans-, and both cis- and trans-regulatory changes between the two species, respectively ().

Classification of transcriptional regulation divergence between ZF and OF.

(A) Heatmaps of gene expression in ZFs and OFs, and allelic expression ratios in F1 hybrids for cis-, trans-, and both cis- and trans-regulated genes in song nuclei HVC and RA (blue–red colored). Comparison between species-different gene expression (A/B) and allelic expression ratios in F1 hybrids (a/b) in heatmaps (dark brown–light yellow colored). “A” and “B” represent RPKM(ZF average) and RPKM(OF average), respectively. “a” and “b” represent Reads (ZF allele) and Reads (OF allele), respectively. (B) Gene numbers classified by cis-, trans-, both cis- and trans-, conserved, and ambiguous regulation in HVC and RA. (C) Scatterplots of expression ratios between ZF and OF (x-axis) and allelic expression ratios in F1 hybrids (y-axis) for genes showing differential expression between species. Blue-, red-, and orange-colored spots: cis-, trans-, both cis- and trans-regulated genes, respectively. Filled spots correspond to species-differentially expressed (SDE) genes. (D) Cis- and trans-effects on the expression of species-differentially regulated genes. The percent of cis-, trans-, both cis- and trans-, conserved, and ambiguous transcriptional regulatory genes in the SDE and non-SDE genes (Fisher’s exact test, ***p < 0.001). Relevant data values are included in . F1, first-generation; OF, owl finch; RA, robust nucleus of the arcopallium; RPKM, reads per kilobase of transcript per million reads mapped; ZF, zebra finch.

In both HVC and RA, trans-alteration was more prevalent than cis-alteration. These results indicated that the expression of 600–800 genes (approximately 10%–15% of the total expressed genes) in the vocal motor song nuclei was modified by altered transcriptional regulation between the two species. Furthermore, a majority of the genes under conserved regulation were highly expressed in both HVC and RA (3,523 genes of 4,489 [78.5%] and 4,782 [73.7%] genes expressed in HVC and RA, respectively). In contrast, most of the cis- and/or trans-regulated genes were not shared between HVC and RA (), showing a brain region–specific transcriptional regulatory alteration. Although this result was obtained based on a cross-replicate comparison of the ASE ratio using eight F1 hybrids, we confirmed this result through an estimation method using the average of ASE of all F1 hybrids [16,40,41], which showed similar rates of cis- versus trans-regulation divergence (see Materials and methods, ).

Cis- and trans-regulatory effects on species-differential expression

We then examined whether the species-differentially expressed (SDE) genes in HVC and RA were affected by the transcriptional regulatory divergence between ZF and OF. Based on the RPKM values of each gene expressed in ZF and OF, 333 and 374 genes showed significantly different expression in HVC and RA, respectively, between the two species (2.9% and 3.3% of the total genes expressed in HVC and RA) (DEseq2 package, p-value corrected by the Benjamini-Hochberg method, p < 0.05; n = 4 each from ZF and OF) (). Totals of 209 and 242 genes of the SDE genes in HVC and RA, respectively, passed the ss-SNPs threshold for calculating the ASE ratio in F1 hybrids. Such SDE genes were significantly enriched with a higher probability of cis-, trans-, and both cis- and trans-regulatory effects compared with those of non-SDE genes, in both HVC and RA (Fisher’s exact test, ***p < 0.001) (). This shows a significant association of transcriptional regulatory changes with SDE genes in the song nuclei.

A predominant effect on cellular molecular function by trans-regulatory divergence

To understand whether transcriptional regulatory divergence has any potential molecular contribution to cellular functions in HVC and RA, we performed GO enrichment analysis using the sets of genes affected by cis-, trans-, and both cis- and trans-regulatory changes. The result showed that more GO categories were enriched for trans-regulated genes compared with the other types of regulatory divergence in both HVC and RA (Fisher’s exact test, p-value corrected by the Benjamini-Hochberg method) (). In particular, we found that GO categories related to neural functions associated with presynapse, chemical synapse transmission, and neuron projection were significantly enriched for RA trans-regulated genes. These results motivated us to focus subsequently on altered trans-regulation in RA.

To predict the potential regulatory mediators driving species differences in the expression of trans-regulated genes in RA, we performed upstream regulatory analyses using Ingenuity Pathway Analysis (IPA) [42,43]. We found that brain-derived neurotrophic factor (BDNF) was the most significant upstream trans-mediator of genes under trans-regulation in RA, which included genes for neural plasticity and dendritic spine formation (glutamate decarboxylase [GAD] 2, NMDA glutamate receptor [GRIN] 2A, neuropeptide Y [NPY], and collapsin response mediator protein [CRMP] 1) (Fisher’s exact test, p = 6.44 × 10−7) (). Amino acid substitution and trans-mediator expression level changes could potentially mediate the trans-regulatory effects to alter downstream gene expression. In line with this prediction of BDNF as a trans-regulatory mediator in RA, we found two amino acid substitutions in BDNF between ZF and OF: Ser45Arg in prodomain and Thr143Met in nerve growth factor (NGF) domain (). Furthermore, BDNF was an SDE gene in HVC between the two species (Student's t test, *p < 0.05) (). In HVC, as an upstream song nucleus connecting to RA, BDNF mRNA is primarily expressed in neurons projecting to RA [44], meaning that HVC could anterogradely secrete BDNF protein to RA via connecting axons as a potential trans-regulation via neural connections. Furthermore, we found differences between species regarding the regulation of the BDNF mRNA expression level in both HVC and RA: OFs had a higher expression level than ZFs at the 3-hour singing condition that induced singing-driven gene expression change, including BDNF () [44,45]. Therefore, in order to uncover the putative trans-regulatory mechanisms of BDNF and to evaluate its potential impacts in generating species-specific songs, we examined how the amino acid substitution and/or expression level of BDNF relates to song structures.

Correlation between individual variations of the species-biased song phenotypes and the BDNF expression level in F1 hybrids

To evaluate the putative trans-regulatory effects mediated by the BDNF amino acid substitution or expression level, we investigated the correlation between song phenotypes and ASE ratio or the expression level of BDNF in F1 hybrids. Considering that neither ZF nor OF are inbred, the interspecies F1 hybrids might present individual variation in the ASE ratio and expression levels of transcribed genes. Consistently, at the transcriptome analysis in F1 hybrids, we realized that F1 hybrids possessed a wide range of individual difference in their ASE ratio and expression level of BDNF mRNA in HVC and RA (), such that each individual F1 hybrid transcribed ZF- and OF-type BDNFs with a unique expression ratio and level. Furthermore, F1 hybrids acquired individually unique songs with a wide range of ZF- and OF-biased features, even though they were reared listening to both ZF and OF songs as models (). We used the same sets of 7 total acoustic and sequential song parameters that showed differences between ZF and OF (5 for acoustic and 2 for sequential parameters) (). We found only one correlation between the ASE ratio of BDNF in RA and the entropy variance of syllables (r = 0.800, p = 0.017, Pearson correlation) (). In contrast, the expression level of BDNF mRNA in HVC had four significant correlations with acoustic and sequential song parameters in F1 hybrids (acoustics: syllable duration [r = −0.862, p = 0.006] and entropy variance [r = −0.822, p = 0.012]; sequence: motif [r = −0.762, p = 0.028] and repetition [r = 0.729, p = 0.040], Pearson correlation) (). These correlational analyses in F1 hybrids point to the BDNF mRNA expression level in HVC (instead of the amino acid substitution) being the most likely RA trans-acting mechanism, which induces anterograde secretion of BDNA protein to RA.

Correlation between individual variation in BDNF expression level and species-biased song structures in F1 hybrids.

(A) Individual variation of BDNF mRNA expression level and ASE ratio between F1 hybrids. (B) Individual variation of learned songs in F1 hybrids that were tutored with ZF and OF songs. Orange solid and blue dotted lines represent the motif and repetitive structure of syllables, respectively. (C) Heatmaps showing the correlation of p-values between the BDNF expression level or ASE ratio and species-biased song phenotypes in F1 hybrids. (D) Correlations between BDNF mRNA expression in HVC and species-biased song structures (syllable duration, entropy variance, motif, and repetition) among F1 hybrid individuals. Relevant data values are included in . ASE, allele-specific expression; BDNF, brain-derived neurotrophic factor; F1, first-generation; OF, owl finch; OZ, F1 hybrid offspring between OF♀ and ZF♂; RA, robust nucleus of the arcopallium; RPKM, reads per kilobase of transcript per million reads mapped; ZF, zebra finch; ZO, F1 hybrid offspring between ZF♀ and OF♂.

Alternation of trans-regulated gene expression and obliteration of learned song features by the pharmacological overactivation of BDNF receptors in the RA

To examine the potential causal links of transcriptional regulation between the BDNF concentration level and the predicted downstream trans-regulated genes, we infused a selective agonist of the BDNF receptor, i.e., tropomyosin receptor kinase B (TrkB), namely 7,8-dihydroxyflavone (7,8-DHF; 10 μg/μL) in vivo, into the RA of adult ZFs by using local retrodialysis () [46]. Transcriptional analysis to compare control (PBS) and 7,8-DHF–infused birds revealed that 570 genes of the 11,655 genes expressed in the RA were differentially identified, with over 4-fold changes between the two groups (DEseq2, p < 0.05) (). Among the differentially expressed 570 genes, 6 of the 21 putative downstream trans-regulated genes of BDNF (shown in ) were found to have an altered expression after the pharmacological activation of the BDNF receptors. This further supports our earlier finding that BDNF could be a potential regulatory mediator of the RA trans-regulated genes.

Obliteration of species specificity of ZF song by BDNF receptor agonist infusion into RA.

(A) Scatterplot indicating RA gene expression in control and 7,8-DHF–infused birds. Dashed lines represent the boundary of the 4-fold expression difference. Darker gray colored dots represent significant differences in expressed genes higher than 4-fold between the control and 7,8-DHF–infused birds. Red colored dots represent downstream trans-regulated genes of BDNF (represented in Fig 5C). (B) Songs before and after infusing BDNF receptor TrkB agonist, 7,8-DHF. Typical examples of songs from control and 7,8-DHF–infused birds. Orange solid lines represent the motif structure of ZF songs. (C) Examples of syllable sequence changes between pre- and post-infusion. Syllable similarity matrices for a pair of songs produced by control and 7,8-DHF–infused birds. (D) Changes in the frequency of motif and repetition in songs at pre- and post-infusion stages (control ZF, n = 3, ZF with 7,8-DHF infusion [7–10 days], n = 5; paired t test, *p < 0.05). Each dot corresponds to individual birds. (E) Examples of syllable acoustic changes between pre- and post-infusion. Scatterplots indicate the distribution of 150 syllables (duration versus mean frequency) from control and 7,8-DHF–infused birds. (F) Changes in syllable acoustics (syllable duration, inter-syllable gap duration, entropy variance, mean FM, and FM variance) of songs at pre- and post-infusion stages (control ZF, n = 3, ZF with 7,8-DFH infusion [7–10 days], n = 5; paired t test, **p < 0.01, *p < 0.05, n.s., not significant). Each dot corresponds to an individual bird. Relevant data values are included in . BDNF, brain-derived neurotrophic factor; FM, frequency modulation; RA, robust nucleus of the arcopallium; RPKM, reads per kilobase of transcript per million reads mapped; TrkB, tropomyosin receptor kinase B; ZF, zebra finch; 7,8-DHF, 7,8-dihydroxyflavone.

Fig 5

Predominant effect on cellular molecular function by trans-regulatory divergence.

(A) GO enrichment analysis of the cis-, trans-, and both cis- and trans-regulated genes in HVC and RA. Size of points represents the number of genes assigned to each GO term. Red lines represent the p-value for significant enrichment (Fisher’s exact test adjusted by the Benjamini-Hochberg method, p < 0.05). (B) Top 7 candidate upstream mediators for trans-regulated genes in RA. (C) Gene–gene connections for BDNF downstream genes. Pink-colored genes are trans-regulated genes in RA. Solid and dotted lines represent directed and undirected regulation, respectively, between connected genes. (D) Comparison of BDNF amino acid sequences between ZF and OF. (E) BDNF mRNA expression level in HVC, RA, and whole brain between ZF and OF at the silent condition based on RNA-seq data. (F) BDNF mRNA expression in the HVC, RA, and the surrounding areas (caudal nidopallium [cN] and archopallium [A], respectively) of ZF and OF at the 3-hour undirected singing condition (n = 4 each). White signals: BDNF mRNA. Scale bars, 1 mm (in left panes) and 200 μm (in right panel). Relevant data values are included in . a.a., amino acid; BDNF, brain-derived neurotrophic factor; GO, Gene Ontology; NGF, nerve growth factor; OF, owl finch; RA, robust nucleus of the arcopallium; RNA-seq, RNA sequencing; RPKM, reads per kilobase of transcript per million reads mapped; ZF, zebra finch.

We also found that song changes following 7,8-DHF infusion, with a lower syllable transition consistency during the early stage (approximately 5 days after drug infusion). In addition, following continuous infusion for up to 2 weeks, adult structured songs gradually became more degraded, leading to the loss of learned song features in adult ZFs (). Although a few of the acoustic parameters (syllable duration and mean FM) maintained the original traits, syllable sequence (i.e., motif and repetitive indexes) and other acoustic parameters (i.e., inter-syllable gap duration, entropy variance, and FM variance) were drastically changed by the infusion of 7,8-DHF (), thus indicating that a precise amount of BDNF contributes to the maintenance of the learned song structures of ZF.

Discussion

Previous studies have demonstrated monogenic effects on adaptive behavioral phenotypes [7,47-49]. In contrast, the genetic basis of polygenic adaptations has been more challenging to pinpoint. Therefore, elucidating various SDE genes and the transcriptional regulatory divergences could be a promising step towards a better understanding of the contribution of multiple genes to the evolution of behaviors. For these, we examined the distribution of cis- and trans-regulatory divergences underlying the differences in gene expression in specific brain regions associated with the production of learned vocalizations between two closely related songbird species.

A number of studies that used entire organ tissues/body showed that there are more significant changes in cis- than trans-regulation between interspecies/lines of fruit flies [14,50], wasps [51], birds [52], and mouse [37]. In contrast, our study revealed that trans-regulatory changes were more prevalent than cis- in determining gene expression differences in song nuclei between two closely related species (). In addition, biological processes associated with neural functions were more enriched for genes showing trans-regulatory divergence in HVC and RA (). This difference in the effects of cis- or trans-regulations on transcriptional divergence could be caused by different methods of estimation using ASE ratio in the F1 hybrids. However, even when we used an estimation method using the average ASE of F1 hybrids, which has the potential to underestimate trans-regulation [40], we obtained a similar result showing that transcriptional regulatory divergence has occurred primarily in trans-regulation. To examine whether the trans-biased regulatory divergence is specific to song nuclei or not, it will be necessary to perform a similar analysis using samples from non-song nuclei or whole brain in songbirds. Furthermore, we set the cutoff with ≥5 reads at each ss-SNP position and median RPKM ≥ 10 to extract genes that were used for the calculation of the ASE ratio. This cutoff threshold is a stricter setting than other published studies [50,51]. Therefore, if we set a lower threshold to extract the ASE genes, the estimated gene number regulated by different transcriptional divergence would be increased.

In principle, two alleles in the cells of an F1 hybrid are regulated in the same trans-regulatory environment. Therefore, differences in expression of two alleles in F1 hybrids should reflect cis-regulatory divergence between the two parental genomes, generating a consistent ASE ratio among F1 hybrids. However, at a considerable number of genes in the song nuclei, we observed a large variance in ASE ratio among F1 hybrids, which we defined as “ambiguously” regulated genes (). Approximately 10% of the total expressed genes in HVC and RA were categorized as “ambiguous” (). Ambiguous regulation could result from intraspecies genomic variation. Indeed, the experimental ZF and OF have not been genetically selected animals. We found 742,302 and 414,040 polymorphic SNPs in the transcribed sequence from the whole brain of ZFs and OFs, respectively (n = 4 each). Therefore, the individual variability in the ASE ratio between F1 hybrids may be caused by intraspecies polymorphisms, which could in turn be additional trans- and cis-regulatory variants underlying individual difference in gene expression in song nuclei.

We found that BDNF is one of several potential upstream mediators for trans-regulated genes in RA (). BDNF transcription, secretion, and actions are directly regulated by neural activity. Secreted BDNF mediates multiple activity-dependent processes, including neuronal differentiation/growth, synapse formation, and plasticity during developmental and adult stages [53-56]. In the song system of songbirds, singing behavior induces BDNF mRNA expression in song nuclei including HVC, suggesting that neural activity-dependent signaling of BDNF regulates neuronal maturation [44,45,57,58]. We had reported that ZFs prevented from singing during the song learning period possess immature dendritic spine density in RA neurons and produced highly unstable song lacking species-specific features when allowed to sing freely, even at the adult stage [59]. Although transient BDNF up-regulation in HVC enhances song learning during the critical period [60], a short-term local injection of BDNF into RA of adult ZFs changed crystallized songs to juvenile-like plastic songs with sequence variability; these changes correlated with an increase in HVC axonal boutons in RA [61]. We further confirmed that the continuous and local infusion of BDNF receptor agonist 7,8-DHF into RA of adult ZFs induced severe song degradation, eliminating both learned acoustic and sequence features. Therefore, we suggest that BDNF mediates the precise synaptic connections and strength of connections allowing HVC to activate populations of RA neurons at specific and precise time points during song rendition. Species differences in BDNF expression in the song nuclei could shape the anatomical and molecular bases for the generation of species-specific learned song structure via the activity-dependent trans-regulation of the downstream genes. Using F1 hybrids, we found that ASE ratios of BDNF were more weakly associated with species-biased song phenotypes compared with BDNF mRNA expression levels. However, this does not rule out a potential trans-regulatory effect mediated by BDNF amino acid substitution on species-specific behaviors. In human, the BDNF polymorphism (Val66Met; rs6265) affects intracellular trafficking and reduces activity-dependent secretion of mature BDNF [62]. The BDNF polymorphism causes altered dendritic spine density, memory formation, and extinction [63,64]. Future application of genome editing technologies would be a powerful tool to elucidate the in vivo contribution of BDNF polymorphisms to species-specific behaviors.

We investigated the divergence between the ZF and OF in terms of gene transcription for the generation of species-specific learned songs (but not for the learning process). Thus, we performed a series of experiments including song comparative analysis, comprehensive RNA sequencing (RNA-seq), and BDNF agonist infusion by using adult birds after the critical period for song learning. However, it is crucial to consider the potential effects of BDNF on the development of neural circuits for species-specific song learning and production during the embryonic and early post-hatching periods. Although we observed that the pharmacological overactivation of BDNF receptors drastically affected song change and led to the loss of learned song structures at both syllable acoustic and sequence levels, we cannot tell whether such song degradation is induced by any species-specific deficiency. Given the wide variety of BDNF cellular functions, the pharmacological experiment was limited by the selective modification of signaling machinery for species-specific song generation. We found that not only predicted downstream trans-regulated genes but also over 550 genes had altered expression levels in the RA, as assessed by comparing control and 7,8-DHF–infused birds. Therefore, future research with more refined experiments for targeted multiple genes, manipulation timing, and cell types will be crucial.

In this study, we investigated the regulatory drivers of species divergence in gene expression in song nuclei in the vocal motor circuit in adults. We suggest these regulatory differences between species could explain a genetic molecular mechanism for the generation and maintenance of the species specificity of learned songs. The anterior forebrain pathway (AFP) is a cortico-basal ganglia-thalamocortical loop, which is a specific pathway for song learning during development and for vocal plasticity maintenance later in life [65-69]. For sensorimotor coordination, AFP generates instructive biased variability and conveys this to the premotor song nuclei RA as a reinforcement signal [46,70]. Currently, we cannot make direct causal links between AFP function and the acquisition of species-specific song patterns. However, lesion of the basal ganglia nucleus, Area X of the striatum (Area X), in the AFP at an early critical period was shown to disrupt motif structure, a sequential trait commonly observed in ZF songs [69]. Furthermore, the expression of transcription factors such as Forkhead box protein P2 (FoxP2) and androgen receptors in Area X shows species-specific patterns [25,26]. These transcription factors could be potential regulators for further species difference due to their regulatory effects on downstream genes, thereby generating species-biased vocal plasticity, which in turn promotes species-specific song learning. Therefore, studying species differences in gene expression in the song nuclei of the AFP through the critical period of song learning would provide vital insight into how species-specific patterns of gene expression underlie species-specific songs.

In conclusion, our results suggest a neurogenetic association between brain region–specific transcriptional divergence and species-specific learned behaviors. Most complex motor skills, such as birdsong and human speech, are acquired through learning and constrained in a species-specific manner [35,71]. Using similar techniques to those developed in the present study on other interspecific hybrids could give additional insights into the existence of more conserved or unique cis-/trans-regulatory changes underlying the evolution of species-specific learned behaviors.

Materials and methods

Ethics statement

All experiments were conducted under the guidelines and approval of the Committee on Animal Experiments of Hokkaido University (Approved No. 18–0053). These guidelines are based on the national regulations for animal welfare in Japan (Law for the Humane Treatment and Management of Animals with partial amendment No. 105, 2011). For brain sampling, the birds were humanely killed by decapitation after injection of an overdose of pentobarbital.

Animals and song tutoring

ZFs (T. guttata) and OFs (T. bichenovii) were obtained from our breeding colony at Hokkaido University and local breeders. Reciprocal F1 hybrids were bred by pairing ZF and OF at our breeding colony. All birds were maintained with food and water available ad libitum under a 13:11-hour light/dark cycle. For song cross-tutoring experiments, ZF chicks were raised by both parents in breeding cages until 10–15 phd, and then the father was removed by 15–25 phd from the cage to prevent male juveniles from listening their father’s song. OF chicks were hand-raised after hatching until they could feed themselves (approximately phd 30–40). After fledging, juveniles were subsequently housed in individual isolation boxes and then individually housed in a sound-attenuating box containing a mirror to reduce social isolation. Cross-species’ tutor songs were played 7 times each in the morning and afternoon at 55–75 decibels from a speaker (SRS-M30, SONY, Tokyo, Japan) passively controlled by Sound Analysis Pro. Similarly, F1 hybrids were song tutored by passively and randomly playing a set of ZF and OF songs with an interval duration of 300–500 ms as the song model.

Song recording and analysis

Songs were recorded using a unidirectional microphone (SM57, Shure) connected to a computer with Sound Analysis Pro (SAP v1.04). For analysis of the acoustic features of songs, 500 syllables were randomly selected from ZF and OF songs (n = 6 birds each). To characterize the syllable that differed between ZF and OF, a total of 10 acoustic features were measured: syllable duration, inter-syllable gap duration, mean pitch, pitch goodness, Wiener entropy, entropy variance, mean AM, AM variance, mean FM, and FM variance [31]. Statistical analysis was performed on these acoustic features between ZF and OF by one-way ANOVA. For the analysis of the sequence feature of songs (motif and repetition rates in a song), a syllable similarity matrix (SSM) analysis was performed following a previously reported method [30] (). This method calculates the contiguous syllables transition frequency of “paired (motif)” and “repetitive” syllables transitions in the songs. To test song structure changes by pharmacological manipulation, we measured both the syllable acoustic and sequential parameters of 150 syllables at pre- and post-time points (7–10 days) after drug infusion. Six acoustic syllable parameters (syllable duration, inter-syllable gap duration, entropy variance, AM variance, mean FM, and FM variance) and sequence features (motif and repetition transition rates) were used for the PCA and 2D view, and this was performed using the prcomp and rgl packages in R, respectively.

Brain tissue sampling and RNA extraction for RNA-seq

For sampling of whole brain tissues, adult male birds were isolated in a soundproof chamber for at least one day before humanely killing (ZF: n = 4, OF: n = 4, adult [>200 phd]). Birds were killed under silent and dark condition in the morning before the lights were turned on. The pallium and pallidum regions were rapidly dissected, frozen on dry ice, and stored at −80°C until RNA extraction. Total RNA was isolated using TRIzol Reagent according to the manufacture’s protocol (Invitrogen) and was treated with RNase-free DNase.

For sampling of HVC and RA tissues by laser capture microdissection (LCM), adult ZF, OF, ZO, and OZ F1 hybrids (n = 4 birds each, >130 phd) and control and 7, 8-DHF–infused ZFs (n = 3 each, >130 phd) were isolated in sound-attenuation boxes and killed under silent and dark condition. Brains were embedded in OCT compound (Sakura Fine Technical) and stored at −80°C until use. Brain sections were cut at a 14-μm thickness in the sagittal plane and mounted onto glass slides with a handmade membrane system for laser microdissection. We confirmed the presence and boundaries of HVC and RA using Nissl staining (LCM Staining kit; Ambion). HVC and RA were microdissected using a laser capture microscope ArcturusXT (Arcturus Bioscience) with the following parameter settings: spot diameter, 100 μm; laser power, 80 mW; and laser duration, 80 ms [72]. The captured tissues were dissolved into RLT buffer (Qiagen) with β-mercaptoethanol, treated with DNase in the column to avoid contamination of genomic DNA, and then stored at −80°C until RNA extraction.

RNA-seq library construction and sequencing

RNA integrity number (RIN) and concentration were measured with Bioanalyzer 2100 (Agilent Technologies) to guarantee the quality of RNA. For RNA-seq of HVC and RA, we performed first-strand cDNA amplification using total RNA (1–2 ng) from HVC and RA under a PCR amplification condition of 14 cycles at 98°C for 10 seconds, 65°C for 15 seconds, and 68°C for 5 minutes, following the Quartz-amplification method [73]. Amplified cDNAs were purified using a PCR purification column (MiniElute PCR Purification Kit; Qiagen) and the concentration was measured using Bioanalyzer 2100 (Agilent Technologies). Non-amplified first-strand cDNAs synthesized using total RNA from the whole brain (telencephalon) and amplified cDNAs using total RNA from HVC and RA tissue were used to construct poly(A) selected paired-end sequencing libraries (TruSeq DNA Sample Prep Kits, Illumina). All libraries were sequenced using the Illumina Hiseq2500 platform for 100-bp paired-end sequencing.

For each telencephalon brain sample, 33.5–47.0 M RNA-seq reads were output from the Illumina Hiseq 2500. Sequencing reads were mapped onto the ZF reference genome obtained from Ensembl (Taeniopygia_guttata taeGut3.2.4.dna.fa) with the Tophat2 program and assembled to predicted transcripts with the Cufflinks program. Through comparison with the previous annotation file using the cuffcompare program, 12,156 transcripts were identified as predicted RNA transcripts expressed in the ZF telencephalon. All RNA-seq data were deposited in the DDBJ Sequence Read Archive (submission numbers DRA005548, DRA002970, and DRA008696).

Identification of ss-SNPs

Adapter sequences of raw data from ZF and OF whole brain NGS results were removed by Trimmonatic. Clean reads from ZF and OF whole brain were mapped to a ZF reference genome obtained from Ensembl (Taeniopygia_guttata.taeGut3.2.4) by TopHat2 to reconstruct pseudo ZF and OF genomes. Mapped reads with longer gaps (>3,000 bp) were removed in the subsequent analysis. ss-SNPs and insert and deletions (indels) between ZF and OF were identified from the mapping result of the whole brain reads. The positions of ss-SNPs and indels were used to reconstruct pseudo genomes of ZF and OF. ss-SNPs were defined as follows: the base variants were same in all individuals of a species, but different from the base found in all individuals of another species. SNPs in individuals of the same species (intraspecies SNPs) were maintained as the same base for both ZF and OF reconstructed genome sequences. MUMmer software was used to identify ss-SNPs using the reconstructed ZF and OF genomes.

Read mapping and quantification of gene expression level

Low-quality reads and adaptor sequences were removed from all HVC and RA raw reads using the Filter FASTQ pipeline (https://cell-innovation.nig.ac.jp) and Flexbar software. Clear reads were mapped to reconstructed ZF genome by TopHat2. Transcript levels were quantified as RPKM value. Cufflinks was used to evaluate the expression levels of each gene by calculating the RPKM of HVC and RA samples of ZF and OF using the improved genome annotation Gene Transfer Format (GTF) file [59]. Based on the RPKM of individual birds (n = 4 each from ZF and OF; n = 3 each from 7,8-DHF–infused and control ZFs), the expression differences of each gene were identified between ZF and OF and between 7,8-DHF and control ZFs as differently expressed genes using the R package DEseq2 (adjusted p-value < 0.05, the Benjamini-Hochberg procedure).

Allelic expression ratio in F1 hybrids

To distinguish reads of the two alleles in F1 hybrids, the mapping results of HVC and RA of F1 hybrids were used following SNPsplit’s instruction. First, an N-marked genome sequence was constructed by replacing “N” at the ss-SNP position in the ZF pseudo genome. RNA-seq reads of HVC and RA of ZF, OF, and F1 hybrids were mapped to the N-marked genome by TopHat2. ss-SNPs were identified as SNP sites having more than 98% of total reads that were different between ZF and OF alleles. The identified ss-SNPs were reverified by reads from HVC and RA of ZF and OF, to enhance the reliability of ss-SNPs. The mapped HVC and RA reads of F1 hybrids were then separated into ZF or OF allele transcripts based on the ss-SNP information, and the number of reads was counted at each ss-SNP position by SAMtools.

The following thresholds were set for calculating the allelic expression ratio of each gene expressed in HVC and RA of F1 hybrids: (i) existence of at least one ss-SNP, (ii) more than 5 reads at each ss-SNP site, and (iii) a median RPKM of at least 10 for all 16 individuals (including ZF, n = 4; OF, n = 4; ZO, n = 4; OZ, n = 4). The allelic expression ratio was quantified using the d-score [74]: d-scores of 0 reflect equal expression between the two alleles, whereas d-scores of −0.5 and 0.5 reflect exclusive transcription from OF or ZF alleles, respectively.

Identification of cis- and/or trans-regulatory divergence

The potential of genomic imprinting in F1 hybrids was tested using Spearman’s rank correlation of gene allelic expression ratio between ZO (n = 4) and OZ (n = 4). The difference in the allelic expression ratio of each gene was compared between ZO and OZ hybrids using one-way ANOVA (ZO, n = 4; OZ, n = 4; adjusted p-value by the Benjamini-Hochberg method).

Cis- and/or trans-regulatory divergences were evaluated using a previously reported method [16]. The gene expression ratio between parental species was calculated with the formula X = log2(A/B), where “X” is the gene expression ratio between parental species; “A” and “B” are the average RPKM for ZFs and OFs, respectively (n = 4 each from ZF and OF). The allelic expression ratio of F1 hybrids was calculated as Y = log2(a/b), where “Y” is the allelic expression ratio between two alleles; “a” and “b” are the read counts of ZF and OF alleles in F1 individuals, respectively. Cis- and trans-effects on gene expression divergence were estimated by the scheme described in . In brief, the regulation mechanism of gene expression between ZF and OF was (1) a cis-regulatory difference if X = Y and Y ≠ 0; (2) a trans-regulatory difference if X ≠ Y and Y = 0; (3) both cis- and trans-regulatory differences if X ≠ Y and Y ≠ 0; (4) no cis- and trans-regulatory differences (i.e., conserved) if X = Y and Y = 0. The Student's t test was used to determine the difference between the gene expression ratio in parental species and the allelic expression ratio in F1 hybrids. The SGoF program was employed to correct p-values for multiple testing (adjusted p ≤ 0.05). The previous standard method for estimating regulatory divergence can lead to a negative correlation as an artifact when cis-estimates have any errors [40, 41]. To avoid this bias, first we randomly selected four individual F1 hybrids as a group to estimate cis-effects using their average ASE ratio while the remaining four F1 individual hybrids were used to compare the expression ratio between ZF and OF. For each gene, a total of 70 combinations were constructed by random selection of four of eight F1 hybrid birds (n = 4 each from ZO and OZ). Thus, cis- and/or trans-regulatory identification was done for each gene for each pair of 70 total combinations. During this cross-replicate comparison, some genes were categorized as different transcriptional regulations due to a large variance in ASE ratios among F1 individuals. Therefore, we finally determined which transcriptional divergence made the main regulatory effect on each gene by two steps of statistics following (i) calculation of the difference between four categories (cis-, trans-, both cis- and trans-, and conserved) using the Chi-squared test (with adjusted p-value by FDR < 0.05) and (ii) a comparison of the difference between the first- and second-strongest regulatory effects using a Fisher’s exact-test (adjusted p-value by FDR < 0.05). If genes did not show significance at both tests, such genes were defined as “ambiguous regulatory genes” ().

In addition, we performed analysis of cis- and/or trans-regulatory divergence using a standard method [37,51] and compared these results with those from the above method. The difference of the standard method is that the allelic expression ratios of all eight F1 hybrids (ZO = 4, OZ = 4) were used to estimate cis- and trans-regulatory effects. In brief, the parental expression ratio value X and the allelic expression ratio in F1 hybrid value Y were calculated similarly to our new method. The average values Y of eight F1 hybrid individuals were compared with values X and 0, respectively, to estimate cis- and trans-effects by the scheme described in Fig 1D (Student's t test). The SGoF program was employed to perform multiple testing correction (adjusted p-value ≤ 0.05) ().

Fig 1

Cis- and/or trans-regulatory changes during species differentiation.

(A) During evolution, cis- and/or trans-regulatory elements change gene expression levels between closely related species. (B) Schematic showing selected song-control regions and connections in the songbird brain. The posterior motor pathway and the anterior cortico-basal ganglia-thalamic circuit (anterior forebrain pathway [AFP]) are represented as red and gray lines, respectively. (C) Genome composition of reciprocal F1 hybrids between zebra finch (ZF) and owl finch (OF). ZO represents F1 hybrid offspring between ZF♀ and OF♂. OZ hybrids are the opposite. Male F1 hybrids share identical sets of auto- and sex chromosomes. (D) Classification of species differences in cis- and/or trans-regulations based on the comparison of the relative gene expression ratio between parental species and the allelic expression ratio in their F1 hybrids. For each gene, “A” and “B” represent gene expression levels in ZF and OF, respectively. “a” and “b” represent gene expression levels from ZF and OF alleles, respectively, in F1 hybrids. “A/B” and “a/b” are the expression ratio between parental species and the allelic expression ratio in F1 hybrids, respectively. Area X, Area X of the striatum; DLM, dorsal lateral nucleus of the medial thalamus; F1, first-generation; HVC, used as a proper name; LMAN, lateral magnocellular nucleus of the anterior nidopallium; nXIIts, tracheosyringeal part of the hypoglossal nucleus; RA, the robust nucleus of the arcopallium.

Functional analysis of cis- and/or trans-regulated genes

The functions of genes with cis-, trans-, and cis- and trans-regulatory divergences between ZF and OF in HVC and RA were annotated by GO analysis (DAVID Bioinformatics Resources 6.8; https://david.ncifcrf.gov). GO enrichment analysis was performed for each gene group using Fisher’s exact tests (p-value was adjusted by the Benjamini-Hochberg method). As trans-regulated genes in RA were enriched for the most GO terms, an upstream regulatory analysis was performed for RA trans-regulated genes using IPA software.

In situ hybridization

BDNF cDNA fragments used for the synthesis of in situ hybridization probes were cloned from a whole-brain cDNA mixture of a male ZF. Total RNA was transcribed to cDNA using Superscript Reverse Transcriptase (Invitrogen) with oligo dT primers. The cDNAs were amplified by PCR using oligo DNA primers directed to the open reading frame region from the NCBI cDNA database. PCR products were ligated into the pGEM-T Easy plasmid (Promega). The cloned sequences were searched using NCBI BLAST/BLASTX to compare with homologous genes to other species and genome loci identified using BLAT of the UCSC Genome Browser.

Adult male ZFs (n = 4) and OFs (n = 4) were used. Birds were individually housed in sound-attenuating boxes overnight. On the following morning, singing behavior (undirected singing) was recorded for 3 hours after the lights were turned on. After each singing behavior observation session, the birds were euthanized by decapitation. Brains were embedded in OCT compound (Sakura Fine Technical) and stored at −80°C until use. Frozen sections (12-μm thick) were cut in the sagittal plane. Brain sections for a given experiment were simultaneously fixed in 3% paraformaldehyde/1× PBS (pH 7.4), washed in 1× PBS, acetylated, dehydrated in an ascending ethanol series, air-dried, and processed for in situ hybridization with antisense 35S-UTP–labeled riboprobes of genes. To generate the riboprobes, gene inserts in the pGEM-T Easy vector were PCR amplified with plasmid M13 forward and reverse primers and then gel purified. The amplified DNA fragments and SP6 or T7 RNA polymerase were used to transcribe the antisense 35S-riboprobes. A total of 1 × 106 cpm of the 35S-probe was added to a hybridization solution (50% formamide, 10% dextran, 1× Denhardt’s solution, 12 mM EDTA [pH 8.0], 10 mM Tris-HCl [pH 8.0], 300 mM NaCl, 0.5 mg/mL yeast tRNA, and 10 mM dithiothreitol). Hybridization was performed at 65°C for 12–14 hours. The slides were washed in 2× SSPE and 0.1% β-mercaptoethanol at room temperature for 1 hour; 2× SSPE, 50% formamide, and 0.1% β-mercaptoethanol at 65°C for 1 hour; and 0.1× SSPE twice at 65°C for 30 minutes each. Slides were dehydrated in an ascending ethanol series and exposed to X-ray film (Biomax MR, Kodak) for 1–14 days. We carefully attended in order not to overexpose X-ray films to S35-riboprobe hybridized brain sections. The slides were then dipped in an autoradiographic emulsion (NTB2, Kodak), incubated for 1–8 weeks, and processed with D-19 developer (Kodak) and fixer (Kodak). For quantification of mRNA signal, exposed X-ray films of brain images were digitally scanned under a microscope (Leica, Z16 APO) connected to a CCD camera (Leica, DFC490) with Application Suite V3 imaging software (Leica), as previously described [45, 72, 75, 76]. To minimize handling bias for signal detection among experimental groups, we performed in situ hybridization using multiple brain sections at once for each probe and exposed S35-riboprobe hybridized brain sections on the same sheet of X-ray films. The same light settings were used for all images. Photoshop (Adobe Systems) was used to measure the mean pixel intensities in the brain areas of interest from sections after conversion to 256 grayscale images.

Pharmacological manipulation

Custom microdialysis probes were built using a microdialysis membrane (SpectralPor, in vivo microdialysis hollow fiber, O.D. = 216 μm; total weight, <0.035 g) attached to a drug reservoir, based on a previously described method [46]. Probes were bilaterally implanted at positions adjacent to RA using stereotaxic coordinates. Before setting the probe, spontaneous neural activity was measured to verify the location of RA. Microdialysis probes were carefully set outside the RA to avoid physically damaging the RA, because damage to the RA could induce song changes. Following surgery, the reservoir was filled every morning with saline until the bird began to sing consistently and its phonological and syntactical features were confirmed not to be damaged by implantation of probe. To ensure the position of microdialysis probes, tetrodotoxin (TTX; 6–12 μM) was infused into the RA in a hemisphere, and a hemi-RA inactivation-induced song change was confirmed. Saline (n = 3 birds) or 7,8-DHF (10 μg/μL in 0.9% NaCl, pH 7.4–7.6, Santa Cruz; n = 5 birds) was then continuously infused during daytime via the injection of approximately 2.5 uL of solution into the outer reservoir of the microdialysis probes 3 to 4 times daily. The manipulated birds were allowed to move freely in a sound-attenuation chamber, and the song of each individual was recorded over 10 days after initiation of drug infusion. The remaining drug volume and infusion speed were checked by using a transparent polyimide tube as the outer reservoir of the microdialysis probes. Probe positioning was evaluated postmortem by histological staining of tissue sections.

ZF and OF species-specific song features.

(A) (Upper panels) SSM analysis for the detection of syllable sequential transition patterns. The SSM comprises two steps: First, a correlation matrix including the syllable similarity scores was prepared using the round-robin comparison of all syllables in two songs to maintain the sequential order of the syllables in the songs. These similarity scores in the matrix were binarized at a threshold at 0.595. Second, the occurrence rate of two patterns of binarized “2 row × 2 column” cells in the SSM was calculated as a percentage of the paired (motif) and repetitive-syllable transition types (see the Materials and methods). (Lower panel) Test examples of the SSM method using artificial song models mimicking the songs with motif and repetitive sequences. (B) The similar distribution range of syllable acoustic traits between ZF and OF. Violin plots of the distribution of syllable duration, inter-syllable gap duration, entropy variance, AM variance, mean FM, and FM variance from ZF and OF that were reared with conspecific song tutoring (total 3,000 syllables from n = 6 birds each and 500 syllables/bird). (C) PCA of the song features of ZFs and OFs reared under conspecific and cross-species song tutoring conditions (“Con”: n = 6 each from conspecific song tutored ZF and OF; “Cross”: n = 4 and 3 from cross-species song tutored ZF and OF, respectively). Relevant data values are included in for panels B and C. AM, amplitude modulation; FM, frequency modulation; OF, owl finch; PCA, principal component analysis; SSM, syllable similarity matrix; ZF, zebra finch.

(TIF)

Click here for additional data file.

Experimental flowchart for the calculation of species-differently expressed genes and characterization of transcriptional regulatory divergence.

(TIF)

Click here for additional data file.

Species differences in gene expression in HVC between ZF and OF.

(Left panels) Expression levels of PRKAA1, NR2E1, and CACNA1E in song nucleus HVC of ZFs and OFs. Gray-colored boxes represent the position of exons for each gene. Dark blue peaks below exons represent read density. (Right panels) Gene expression levels in ZF and OF. Each dot represents RPKM value for the individual. Mean ± SEM (n = 4 birds each, one-way ANOVA, *p < 0.05; n.s., not significant). Relevant data values are included in . CACNA1E, calcium voltage-gated channel subunit alpha 1E; NR2E1, nuclear receptor subfamily 2 group E member 1; OF, owl finch; PRKAA1, protein kinase AMP-activated catalytic subunit alpha 1; RPKM, reads per kilobase of transcript per million reads mapped; ZF, zebra finch.

(TIF)

Click here for additional data file.

No genomic imprinting genes in reciprocal F1 hybrids of ZF and OF.

Scatterplots of allelic expression ratios of 5,849 and 6,328 genes in HVC and RA, respectively, of OZ and ZO hybrids (Spearman correlation coefficient). Relevant data values are included in . F1, first-generation; OF, owl finch; OZ, F1 hybrid offspring between OF♀ and ZF♂; RA, robust nucleus of the arcopallium; ZF, zebra finch; ZO, F1 hybrid offspring between ZF♀ and OF♂.

(TIF)

Click here for additional data file.

Cis-, trans-, both cis- and trans-, and conserved regulation in HVC and RA estimated by a method using the average of ASE of all F1 hybrids.

Relevant data values are included in . ASE, allele-specific expression; F1, first-generation; RA, robust nucleus of the arcopallium.

(TIF)

Click here for additional data file.

SDE genes in HVC and RA.

(A) SDE genes in HVC and RA. Orange- and blue-colored spots represent significantly higher expression in ZF or OF, respectively (DEseq2 corrected with the Benjamini-Hochberg method, p < 0.05). (B) Venn diagram representing the number of genes in HVC and RA that are differently expressed between ZF or OF. Relevant data values are included in . OF, owl finch; RA, robust nucleus of the arcopallium; SDE, species-differentially expressed; ZF, zebra finch.

(TIF)

Click here for additional data file.

Gene–gene connections driven by the top 7 candidate upstream mediators for trans-regulated genes in RA.

Top 7 candidate upstream mediators, including BDNF, HTT, POU3F1, MAPT, MNKK1, PSEN1, and HDAC4. Trans-regulated genes by BDNF in RA are noted in red. Orange- and green-colored genes are trans-regulated genes that are significantly expressed more highly in RA of ZF or OF, respectively. Relevant data values are included in BDNF, brain-derived neurotrophic factor; HDAC4, histone deacetylase 4; HTT, huntingtin; MAPT, microtubule-associated protein tau; MNKK1, MAP kinase-interacting serine/threonine protein kinase 1; OF, owl finch; POU3F1, POU class 3 homeobox 1; PSEN1, presenilin 1; RA, robust nucleus of the arcopallium; ZF, zebra finch.

(TIF)

Click here for additional data file.

Untethered microdialysis for pharmacological manipulation of BDNF receptors in RA.

(Left) Photograph of homemade microdialysis probe. (Right) A ZF with microdialysis probes bilaterally implanted in RA. BDNF, brain-derived neurotrophic factor; RA, robust nucleus of the arcopallium; ZF, zebra finch.

(TIF)

Click here for additional data file.

Underlying data for Fig 2B and 2C and S1 Fig.

(XLSX)

Click here for additional data file.

Underlying data for Figs 3A and 2B and S3 Fig.

(XLSX)

Click here for additional data file.

Underlying data for Fig 4A–4D and S4–S6 Figs.

(XLSX)

Click here for additional data file.

Underlying data for Fig 5A, 5B, 5E and 5F and S7 Fig.

(XLSX)

Click here for additional data file.

Underlying data for Fig 6A, 6C and 6D.

(XLSX)

Click here for additional data file.

Underlying data for Fig 7A and 7C–7F.

(XLSX)

Click here for additional data file.

23 Jul 2019

Dear Dr Wada,

Thank you for submitting your manuscript entitled "Transcriptional regulatory divergence underpinning species-specific learned vocalization in songbirds" for consideration as a Research Article by PLOS Biology.

Your manuscript has now been evaluated by the PLOS Biology editorial staff, as well as by an academic editor with relevant expertise, and I'm writing to let you know that we would like to send your submission out for external peer review.

However, before we can send your manuscript to reviewers, we need you to complete your submission by providing the metadata that is required for full assessment. To this end, please login to Editorial Manager where you will find the paper in the 'Submissions Needing Revisions' folder on your homepage. Please click 'Revise Submission' from the Action Links and complete all additional questions in the submission questionnaire.

**Important**: Please also see below for further information regarding completing the MDAR reporting checklist. The checklist can be accessed here: https://plos.io/MDARChecklist

Please re-submit your manuscript and the checklist, within two working days, i.e. by Jul 25 2019 11:59PM.

Login to Editorial Manager here: https://www.editorialmanager.com/pbiology

During resubmission, you will be invited to opt-in to posting your pre-review manuscript as a bioRxiv preprint. Visit http://journals.plos.org/plosbiology/s/preprints for full details. If you consent to posting your current manuscript as a preprint, please upload a single Preprint PDF when you re-submit.

Once your full submission is complete, your paper will undergo a series of checks in preparation for peer review. Once your manuscript has passed all checks it will be sent out for review.

Feel free to email us at plosbiology@plos.org if you have any queries relating to your submission.

Kind regards,

Roli Roberts

Roland G Roberts, PhD,

Senior Editor

PLOS Biology

==================

INFORMATION REGARDING THE REPORTING CHECKLIST:

PLOS Biology is pleased to support the "minimum reporting standards in the life sciences" initiative (https://osf.io/preprints/metaarxiv/9sm4x/). This effort brings together a number of leading journals and reproducibility experts to develop minimum expectations for reporting information about Materials (including data and code), Design, Analysis and Reporting (MDAR) in published papers. We believe broad alignment on these standards will be to the benefit of authors, reviewers, journals and the wider research community and will help drive better practise in publishing reproducible research.

We are therefore participating in a community pilot involving a small number of life science journals to test the MDAR checklist. The checklist is intended to help authors, reviewers and editors adopt and implement the minimum reporting framework.

IMPORTANT: We have chosen your manuscript to participate in this trial. The relevant documents can be located here:

MDAR reporting checklist (to be filled in by you): https://plos.io/MDARChecklist

**We strongly encourage you to complete the MDAR reporting checklist and return it to us with your full submission, as described above. We would also be very grateful if you could complete this author survey:

https://forms.gle/seEgCrDtM6GLKFGQA

Additional background information:

Interpreting the MDAR Framework: https://plos.io/MDARFramework

Please note that your completed checklist and survey will be shared with the minimum reporting standards working group. However, the working group will not be provided with access to the manuscript or any other confidential information including author identities, manuscript titles or abstracts. Feedback from this process will be used to consider next steps, which might include revisions to the content of the checklist. Data and materials from this initial trial will be publicly shared in September 2019. Data will only be provided in aggregate form and will not be parsed by individual article or by journal, so as to respect the confidentiality of responses.

Please treat the checklist and elaboration as confidential as public release is planned for September 2019.

We would be grateful for any feedback you may have.

22 Aug 2019

Dear Dr Wada,

Thank you very much for submitting your manuscript "Transcriptional regulatory divergence underpinning species-specific learned vocalization in songbirds" for consideration as a Research Article by PLOS Biology. Your manuscript was evaluated by the PLOS Biology editors as well as by an Academic Editor with relevant expertise and by three independent reviewers.

Based on the reviews, we will probably accept this manuscript for publication, assuming that you remove the BNDF agonist injection experiments as requested by reviewer 2 and modify the manuscript to address all other concerns raised by the reviewers.

We expect to receive your revised manuscript within two weeks. Your revisions should address the specific points made by each reviewer. In addition to the remaining revisions and before we will be able to formally accept your manuscript and consider it "in press", we also need to ensure that your article conforms to our guidelines, one of which is described below under DATA POLICY and is marked with "***IMPORTANT: ". A member of our team will be in touch shortly with a set of requests. As we can't proceed until these requirements are met, your swift response will help prevent delays to publication.

Upon acceptance of your article, your final files will be copyedited and typeset into the final PDF. While you will have an opportunity to review these files as proofs, PLOS will only permit corrections to spelling or significant scientific errors. Therefore, please take this final revision time to assess and make any remaining major changes to your manuscript.

Please note that you may have the opportunity to make the peer review history publicly available. The record will include editor decision letters (with reviews) and your responses to reviewer comments. If eligible, we will contact you to opt in or out.

Please note that an uncorrected proof of your manuscript will be published online ahead of the final version, unless you opted out when submitting your manuscript. If, for any reason, you do not want an earlier version of your manuscript published online, uncheck the box. Should you, your institution's press office or the journal office choose to press release your paper, you will automatically be opted out of early publication. We ask that you notify us as soon as possible if you or your institution is planning to press release the article.

To submit your revision, please go to https://www.editorialmanager.com/pbiology/ and log in as an Author. Click the link labelled 'Submissions Needing Revision' to find your submission record. Your revised submission must include a cover letter, a Response to Reviewers file that provides a detailed response to the reviewers' comments (if applicable), and a track-changes file indicating any changes that you have made to the manuscript.

Please do not hesitate to contact me should you have any questions.

Sincerely,

Di Jiang, PhD

Associate Editor

on behalf of

Roland G Roberts, PhD,

Senior Editor

PLOS Biology

------------------------------------------------------------------------

DATA POLICY:

You may be aware of the PLOS Data Policy, which requires that all data be made available without restriction: http://journals.plos.org/plosbiology/s/data-availability. For more information, please also see this editorial: http://dx.doi.org/10.1371/journal.pbio.1001797

Note that we do not require all raw data. Rather, we ask that all individual quantitative observations that underlie the data summarized in the figures and results of your paper be made available in one of the following forms:

1) Supplementary files (e.g., excel). Please ensure that all data files are uploaded as 'Supporting Information' and are invariably referred to (in the manuscript, figure legends, and the Description field when uploading your files) using the following format verbatim: S1 Data, S2 Data, etc. Multiple panels of a single or even several figures can be included as multiple sheets in one excel file that is saved using exactly the following convention: S1_Data.xlsx (using an underscore).

2) Deposition in a publicly available repository. Please also provide the accession code or a reviewer link so that we may view your data before publication.

***IMPORTANT: Regardless of the method selected, please ensure that you provide the individual numerical values that underlie the summary data displayed in the following figure panels: Figure 2BC, 3AB, 4ACD, 5ABE, 6ACD (this figure might be removed depending on how you meet the request that you must remove the BNDF agonist injection experiments), S1, S3, S4, S6, S7, as they are essential for readers to assess your analysis and to reproduce it. ***IMPORTANT: Please also ensure that figure legends in your manuscript include information on where the underlying data can be found. You may say in every relevant figure legend that, e.g., "Data associated with this figure can be found in the supplemental data file (S1 Data).".

Please ensure that your Data Statement in the submission system accurately describes where your data can be found.

------------------------------------------------------------------------

Reviewer remarks:

Reviewer #1 (Yisi Zhang, signed review):

In this work, Wang et al. studied the difference between two closely related but behaviorally different songbird species, the zebra finch and the owl finch, of their transcriptional regulation in the singing system. In particular, the authors investigated the cis- and trans-regulation within HVC and RA and found that the trans-regulation was more prevalent in both areas. The cis- and trans-regulations were closely related to the species-differential gene expression, especially the trans-regulated genes in RA mediated by BDNF. The latter was differentially expressed in HVC. The widespread BDNF expression levels in HVC among the F1 hybrids were significantly correlated with some song features. Over-activation of the BDNF receptor in zebra finch RA eliminated the species-specific song features and altered the expression of some genes in RA. This is an interesting and well-designed comparative study elucidating the evolutional variation of transcriptional regulation within the song production system using RNA sequencing and pharmacological manipulation. The analysis is comprehensive with sufficient statistical analysis. The arguments will be strengthened if the following issues can be clarified.

1. How different are the genomes of these two species? Since the SNPs are mapped to the ZF reference, does this method naturally introduce bias to the detection of ss-SNPs and the estimation of differential expression? Please discuss.

2. The authors hypothesized that BDNF regulates the trans-regulated genes in RA through axon terminal secretion from HVC to RA. Although the infusion of 7,8-DHF altered song structure, the change looks a lot like a song degradation as opposed to a shift to OF-like song. The hypothesis will be strengthened if there is a difference in the concentration of RA BDNF protein between the two species and/or the F1 hybrids.

3. The authors pointed out that there is no difference in the ASE ratios of the ZO and OZ hybrids and treated them as equivalent. In Figure 6, however, there seems to be a clear division in the BDNF expression levels (at least in RA) between the ZO and OZ. I think this is worth more elaboration or discussion, especially it is interesting that there seems to be more influence from the female side (if I follow the labeling convention of figure 1). I am curious to see if there is a global difference in the song features between ZO, OZ, as well as a comparison with the ZF and OF tutor songs.

There are some minor issues:

-Ln107, “Fig. 2”. Please specify which sub-figures one should look at. Also, it is not very clear what Fig. 2B is showing, and it is very hard for people not in the songbird field to understand the similarity matrices without explanation.

-Why these particular acoustic features were selected but not, for example, entropy, pitch or duration variance?

-Ln173, “Fig. 4B” should be “4A”.

-Ln258, “diurnal singing in HVC neurons” typo.

-Ln314, “are commonly regulated the trans-regulatory environment” typo.

-Ln402, “OF checked” typo.

-Lns901-906, figure #s are not ordered.

Reviewer #2:

Review of “Transcriptional regulatory divergence underpinning species-specific learned vocalization in songbirds”

This is a well-written paper form an investigator that has done important contributions to the biology of bird song. In this paper they do pioneering work as they try to identify the genetic basis that could account for differences in the types of song produced by different species.

The main strength of the paper is that the strategy is novel for songbirds. A similar approach has been used by Hopi Hoekstra for different behaviors in closely related rodent species, but for birds this is clearly the first such study. The experiments are well-performed, and the explanations are clear. My main concern is that similarly to the Hokestra papers, it is not realistic to expect to be able to pinpoint that the differences for complicated behaviors are due to changes in a single gene. Fundamentally, the final interpretation of this type of approach is that these behavioral differences are due to multiple genes interacting in complex ways. Unfortunately, some features of biology cannot be reduced to simple mechanistic models, and this work is a good example of that caveat.

My main criticism about the paper is that I do not believe that the interpretation of the BDNF receptor manipulation is correct. Fundamentally, when they inject the BDNF receptor agonist they see a “simplification” of the song of the zebra finch song – the quality of the song is degraded and it becomes less elaborate. This does not mean that it loses “species-specific acoustic and sequence traits”, it simply means that it becomes degraded. Many other experiments have song that different lesions of the song circuit produce simplified songs that are quite stereotypical. A control for this experiment would be to inject into the birds any other type of drug that affects genes that are NOT present in the genes that they identify as being differentially regulated between zebra and owl finches, and see the behavioral effects. For example, if they inject a drug that messes up a biological process not directly relevant to song production (for example, tubulin polymerization), my guess is that they will see a behavioral change very similar to what they observe with the BDNF agonist (note: they should choose a drug that is not directly toxic to the neurons, but it just messes them up in a reversible manner).

In a more general way, I do not find the logic of this experiment very convincing. The production of song is an extremely complex process that depends on the formation of connection between different brain nuclei throughout development. The manipulation that they do in this paper involves injecting a BDNF receptor agonist in the circuit of a fully formed circuit. The specific connections between the different brain nuclei have formed during the embryonic period and during song learning. However, activating BDNF receptor in an adult animals is just “messing” with up the connections that were already there. It is interesting that the song gets degraded in a more or less predictable manner, but calling this “species-specific features is not justified.

In summary, I think that paper is well done, and is interesting, but I do not agree that the interpretation of the BNDF experiments should be toned down drastically – BDNF receptor activation clearly changes the song, but we cannot tell that it is changing with any type of specificity. Actually, I think that the paper would be stronger without the BNDF agonist injection experiments.

More specific comments:

. Line 209: In particular, we found that GO categories related to neural functions associated with synapse transmission, somatodendritic compartment, and nervous system development were significantly enriched for RA trans-regulated genes.

I found this statement tautological –genes associated with“synapse transmission, somatodendritic compartment, and nervous system development” comprise a huge number of the genes expressed in neurons. Thus, it is expected that simply because of their sheer number, those genes are the ones whose pattern of expression are going to change.

. Line 221 “ this prediction of BDNF as a trans-regulatory mediator in RA, we found two amino acid 
 substitutions in BDNF between ZF and OF: Ser45Arg in prodomain and Thr143Met in NGF domain (Fig. 5D)”.

And Line 247: “In contrast, the expression level of BDNF mRNA in HVC had four 
significant correlations with acoustic and sequential song parameters in F1 hybrids [acoustics: syllable duration (r = −0.862, p =0.006) and entropy vibrance (r = −0.822, p = 
 0.012); sequence: motif (r = −0.762, p = 0.028) and repetition (r = 0.729, p =0.040), 
 Pearson correlation] (Fig. 6C, D). “

I find this argument very confusing. The main type of measurements that they do are related to levels of gene expression. However, if in addition to changes in levels of expression, there are also changes in the sequence of the genes this would make it extremely difficult to separate which of these 2 factors is responsible for the behavioral changes that they observed.

Reviewer #3:

This paper uses F1 species hybrids to investigate how transcriptional regulatory divergence between species relates to the production of species-specific songs. I generally found the question to be interesting and well-addressed. I have a few concerns, primarily about the conclusions reached with the last experiments, and well as some more minor editorial comments. In general, there are a number of places where the manuscript needs more editing for clarity and readability.

My main concern is with the conclusion that over-activation of BDNF eliminates species-specific acoustic and sequence traits. While there are clearly changes to ZF song with over-activation of BDNF, those changes are reminiscent of simplification of song structure and increases in repeats seen with a range of other manipulations where you lose syllable complexity/diversity and increase repetition. For example, the greater stereotypy and syllable degradation you see with lesioning LMAN/DLM in young birds, or changes with deafening or with different acoustic experience during development. While a number of measures of song move in the direction of the Owl Finch song, I’m not convinced that is the same as the song becoming more like OF sing vs. just becoming simpler and noisier. It’s also not clear why changes to BDNF level would make the song more like OF song, unless it was thought that the shared ancestor of OF and ZF had a song more similar to OFs, which isn’t clear based on the songs of other grassfinch species in the clade or something discussed here.

I am also surprised by the weak effect of learning shown in Figure 1. Because 15-25 phd is really late for removing the father, I wonder if this is why there is not a stronger effect of the cross tutoring?

Drug infusion: 10 days is a long time for a permanently implanted microdialysis probes, were there any issues with gliosis? How was it “infused continuously during daytime”? Is it using the reservoir system? If so, Is there any data to demonstrate that there is still drug moving from the reservoir into the probe itself at the end of the day? For example, when was song recorded for the TTX controls? If it is late in the day, that could support that drug is continuing to be infused.

Minor comments

Line 42 analog not homologue

Line 47 eliminated species-specific acoustic and sequence traits is too strong.

Line 54-55 I’m not sure that ecology makes sense here. Not clear on what this is saying

Line 70 song system? Song pathway makes it sound especially linear

Line 80 refs are mostly about Area X

Line 91-92 this sentence is difficult to follow, rephrase

Line 92 define Gene Ontology

Line 102 Is this really a concern, that lab rearing would eliminate species-specific song features? There needs to be better justification for this.

Line 108 I’m not convinced that having somewhat overlapping values for each individual acoustic feature is an indication that there are not some physical constraints of peripheral vocal organs. Also, the sentence is confusing, I think it means to say that because the distribution ranges are overlapping, they can confirm there is not a physical constraint, but that it not how it is written on Line 109.

Line 116-120 It looks like the statistics were to compare within each tutoring condition (normal vs. cross-species). I’m not clear on what was done statistically to say that there are song tutoring effects on most of the song parameters. This would mean whether there is a difference between ZF tutored ZFs and OF tutored ZFs, correct? But none of those statistics are reported. The differences within the cross-tutored groups are convincing, and the main point of the figure, but as written this paragraph makes me wonder about the effects of tutoring. I think a statistical model with tutoring condition, species, and the interaction term would be useful here. It also seems like a lot of statistical tests (8 different song measures) to do with a very small sample size.

Line 166-167 “we identified that over 75 and 10% of examined genes” I don’t understand what this means.

Line 183 what does “between a pair of eight F1 hybrids” mean?

Line 240-243 This statement needs to be more clear. Birds were reared hearing both ZF and OF songs? Just single examples of each? But each bird acquired an individually-unique song?

Line 245 I don’t think it’s “as a result”, delete this phrase.

Line 249 Entropy variance? Not “vibrancy”

Line 258-261 The justification for looking at BDNF needs to be explained more clearly.

Line 314 word missing, regulated ___ in? by?

Line 342 axonal boutons, not buttons

Line 347-349 This sentence needs rephrasing

Line 361-363 Needs a better description of the AFP. The basal ganglia is homologous across vertebrates. The AFP is a specialized cortical-basal ganglia-thalamic loop

Line 369 I’m not sure that “motif structure” is a species specific trait of ZF songs. This seems to be referring to a very specific feature of motif structure, but that’s not clearly defined here.

Line 402 OF chicks not checked.

Line 581 rephrase, was infused to one hemisphere of RA.

Fig 2 Plots in 2B should be on the same scale for con and cross

Fig 6 Labels in legend are wrong relative to the figure (also out of order)

In general, there needs to be more information in the Methods about the rearing of F1 hybrids and how they were tutored.

5 Sep 2019

Submitted filename: PLoS Biology Response to 1st decision v5.2.docx

Click here for additional data file.

18 Sep 2019

Dear Dr Wada,

On behalf of my colleagues and the Academic Editor, Asif A. Ghazanfar, I am pleased to inform you that we will be delighted to publish your Research Article in PLOS Biology.

The files will now enter our production system. You will receive a copyedited version of the manuscript, along with your figures for a final review. You will be given two business days to review and approve the copyedit. Then, within a week, you will receive a PDF proof of your typeset article. You will have two days to review the PDF and make any final corrections. If there is a chance that you'll be unavailable during the copy editing/proof review period, please provide us with contact details of one of the other authors whom you nominate to handle these stages on your behalf. This will ensure that any requested corrections reach the production department in time for publication.

PRESS

We frequently collaborate with press offices. If your institution or institutions have a press office, please notify them about your upcoming paper at this point, to enable them to help maximise its impact. If the press office is planning to promote your findings, we would be grateful if they could coordinate with biologypress@plos.org. If you have not yet opted out of the early version process, we ask that you notify us immediately of any press plans so that we may do so on your behalf.

We also ask that you take this opportunity to read our Embargo Policy regarding the discussion, promotion and media coverage of work that is yet to be published by PLOS. As your manuscript is not yet published, it is bound by the conditions of our Embargo Policy. Please be aware that this policy is in place both to ensure that any press coverage of your article is fully substantiated and to provide a direct link between such coverage and the published work. For full details of our Embargo Policy, please visit http://www.plos.org/about/media-inquiries/embargo-policy/.

Thank you again for submitting your manuscript to PLOS Biology and for your support of Open Access publishing. Please do not hesitate to contact me if I can provide any assistance during the production process.

Kind regards,

Hannah Harwood

Publication Assistant,

PLOS Biology

on behalf of

Roland Roberts,

Senior Editor

PLOS Biology