Five black-capped chickadees and eight American goldfinches were captured at Concord Field Station in Boston, Massachusetts, using mist nets and cage traps. Similarly, three black-capped chickadees and three American goldfinches were captured at a field station near Kent State University in Ravenna, Ohio, using mist nets. The birds were selectively captured at the area of occurrence, based on specific criteria, such as age, sex, or breeding status. We ensured the swift release of unsampled birds into the wild to reduce the duration of captivity and minimize the disruption to flock dynamics. This approach ensured that non-breeding and non-dominant birds in the flock were only sampled for this study.

To ensure consistency, all captured birds were subsequently housed in cages for 24 hours at room temperature, maintaining a 12-hour light photoperiod. During this time, they were provided with unrestricted access to both food and water. This standardized procedure was implemented to normalize stress levels among all the captured birds before tissue harvest.

The birds were humanely euthanized using the rapid cardiac compression method as per approved Kent State University IACUC protocol (see Ethics Statement below). Birds were immediately dissected to obtain the brain tissues and preserved in RNALater. Samples were further stored at -80°C, until downstream processing. Approximately 25 mg of brain tissue was used for the extraction of total RNA, using the QIAGEN RNeasy Tissue protocol. Following quality control procedures, the samples that met the required standards were sent for mRNA sequencing using the Illumina NovaSeq platform.

We used Trinity v2.9.1 (Grabherr et al., 2011; Haas et al., 2013) to build de novo reference transcriptome assembly using RNAseq data from one representative individual of each population. To assess the completeness of the transcriptome assemblies, the raw RNAseq reads were aligned back to the final transcriptomes. using Bowtie2 version 2.5.0 (Langmead et al., 2009). To further evaluate the completeness of the assembled transcriptomes, we estimated the contig Nx statistic that calculated the length of the shortest contig at which all contigs the same length or longer cover at least 50% of the genome assembly. An expression dependent N50 statistic was also calculated by factoring in the gene expression data, which is a better indicator of transcriptome assembly quality than the Nx statistic alone. The quality and completeness of the transcriptome assemblies were further assessed using BUSCO v5.2.2 (Manni et al., 2021). For this analysis, the set of conserved genes in Passeriformes was selected as a reference from the OrthoDB v10 database (Kriventseva et al., 2019).

We appreciate that the sample size of these bird populations varies (black-capped chickadee - five from Boston population and three from Kent population; American goldfinch – eight from Boston population and three from Kent population). Ideally, we would have wanted to balance the same size in these populations, however the fieldwork was conducted during the time of global covid pandemic, and it was challenging to manage logistics in the field to obtain an ideal sample size in these two populations. To address this unequal sampling issue, we utilized a method that normalize RNA sequencing read counts across samples and control for sample size bias using the tool DESeq2 (Love et al., 2014). DESeq2 uses the median of ratios method for normalization, accounting for library size differences and compositional biases in RNA-seq data. By calculating the geometric mean of each gene across all samples and then dividing the counts for each gene in each sample by this mean, DESeq2 aims to normalize the data while accounting for differences in sequencing depth and RNA composition. This normalization procedure helps mitigate biases that can arise due to varying expression levels of different genes across conditions or sample size biases.

We then used Kallisto v0.46.1 (Bray et al., 2016) on the normalized RNAseq data to quantify the abundance of transcripts in each sample, followed by Sleuth v0.30.0 (Pimentel et al., 2017), to conduct statistical analysis to identify differentially expressed genes, using a cutoff q-value (Benjamini-Hochberg multiple testing corrected p-value) < 0.1 (Pimentel et al., 2017). The interactive web interface (Sleuth Live) was further used to generate clustered sample heatmaps, conduct Principal Component Analyses (PCA) and generate volcano plots for displaying differential expression results among groups.

The Trinotate v3.2.2 pipeline (Bryant et al., 2017) was used for the annotation of the transcriptomes (see Supplementary Text for details). PANTHER v18.0 database (Mi el al., 2013), STRING v12.0 (Szklarczyk et al., 2023) and g:Profiler was used for functional analysis of the differentially expressed genes (see Supplementary Text for details).

We mapped the RNAseq reads to high-quality reference genomes of black-capped chickadee (Poecile atricapillus) (Sayers et al., 2022) and the American goldfinch (Spinus tristis) (manuscript in preparation) using GATK-recommended best practices for variant discovery with RNA-Seq data (Brouard et al., 2019). We used a recently developed workflow from GATK (https://github.com/gatk-workflows/gatk4-rnaseq-germline-snps-indels), that is specifically optimized to call SNPs from RNAseq data to address challenges of calling SNPs from RNAseq data, arising from unequal expression of alleles (see Supplementary Text for details).

We further utilized SNPEff v.5.2a (Cingolani et al., 2012) to conduct functional annotation of the filtered SNPs generated above, using genome annotation database from both species. We predicted coding effects such as synonymous or non-synonymous amino acid changes, start codon gains or losses, stop codon gains or losses, or frame shifts.

We conducted a comparison between two sets of genes: one comprising genes that exhibited significant differences in expression levels between the Boston and Kent populations in both species and the other containing genes that shared regions with SNPs displaying significant differences in allele frequencies between these same populations. The goal was to determine if there was any functional overlap between these two sets of genes. To do this, we employed STRING v12.0 to explore potential functional interactions between the proteins encoded by these two gene sets. We generated protein-protein interaction networks that encompassed all interconnected nodes, with each node representing all proteins produced by a single protein-coding locus. Within these networks, the edges were indicative of both the type and strength of evidence for protein-protein associations between genes.

The sampling locations’ longitude and latitude coordinates were utilized in Daymet’s Single Pixel Extraction Tool (Thornton et al., 2020) to retrieve minimum temperature, maximum temperature, and snow water equivalent for specified date ranges (see Supplementary Text for details).

We captured black-capped chickadees and American goldfinches from two locations (Kent and Boston) and collected their brain tissue samples ( Table 1 ).

RNA was obtained from each of these samples ( Supplementary Table 1 ). The RNA integrity number (RIN) score was 9.5 ± 0.20 (on a scale of 1-10) for all samples, which indicated high quality. RNA sequencing yielded an average of 46 million reads per sample ( Supplementary Table 1 ).

De-novo transcriptome assemblies for each species were generated using these RNAseq reads. The black-capped chickadee transcriptome assembly had a total of 144,272 unique transcript contigs, with a Contig N₅₀ of 3,456 bp ( Table 2 ). Similarly, The American goldfinch transcriptome assembly had a total of 180,455 unique transcript contigs, with Contig N₅₀ of 3,417 bp. RNAseq reads alignment rate to the reference transcriptome was ~95% for black-capped chickadee and ~96% for American goldfinch, indicating that a majority of the raw RNA-Seq reads were utilized to create the transcriptome assembly for both species.

BUSCO analysis indicated that black-capped chickadee transcriptome had a total of 80.9% complete BUSCOs, with 30.3% complete and single-copy, 50.6% complete and duplicated, 2.2% fragmented, and 16.9% missing out of 10,844 total BUSCO groups searched ( Supplementary Figure 1 ). The American goldfinch transcriptome had similar statistics (a total of 80.0% complete BUSCOs, with 27.3% complete and single-copy, 52.7% complete and duplicated, 2.3% fragmented, and 17.7% missing out of 10,844 total BUSCO groups searched ( Supplementary Figure 1 ).

We identified 108 differentially expressed genes (DEGs) in brain tissues between black-capped chickadee (BCCH) populations from Boston and Kent ( Figure 3 ). However, in the case of American goldfinch (AMGO) populations from Boston and Kent, we only identified 14 DEGs ( Figure 3 ). Out of 108 DEGs in BCCH, 96 were annotated to a known gene (80 were annotated using the Trinotate annotation pipeline and 16 were annotated manually, Supplementary Table 2 ). Similarly, for AMGO, 11 out of 14 genes were annotated to a known gene ( Supplementary Table 3 ). We further examined the expression of these DEGs across all samples ( Figure 4 ; Supplementary Figure 2 ), which showed consistent differences between all individuals of Boston and Kent populations (either up or downregulated) in both species. The principal component analysis based on the transcript abundance data showed two clusters consistent with sample geographical location, with PC1 and PC2 explaining 96.27% of the total variance in BCCH ( Figure 5 ) and 97.76% of the total variance in AMGO ( Supplementary Figure 3 ).

The PANTHER v18.0 database was used to extract the gene family and protein class information for these DEGs in black-capped chickadees ( Supplementary Table 4 ) and American goldfinches ( Supplementary Table 5 ). We further carried out gene functional enrichment analysis of these DEGs. For the set of black-capped chickadee DEGs upregulated in the Kent population, STRING detected one enriched Gene Ontology (GO) term whereas g:Profiler detected functional enrichment in seven GO terms ( Table 3 ). For the set of black-capped chickadee DEGs upregulated in Boston, STRING did not identify any enrichment in gene functions. g:Profiler, however, detected enrichment for one GO term ( Supplementary Table 6 ). For AMGO, neither STRING nor g:Profiler identified any enrichment in GO terms among the DEGs.

On average 82.43% of black-capped chickadee (BCCH) RNAseq reads and 85.39% of American goldfinch (AMGO) reads across all samples were mapped to their respective reference genomes.

These alignments were further used to identify 1,538,703 sites that were polymorphic in at least one sample in BCCH. However, we only selected 433,210 high-quality biallelic SNPs in subsequent analyses using rigorous variant filtering criteria, as outlined in the methods. Using similar methods, we identified 535,179 high-quality biallelic SNPs that were polymorphic in AMGO.

For black-capped chickadees, we identified 75,780 non-synonymous SNPs and 247,643 synonymous SNPs (N/S ~ 0.3). Out of 17,929 genes, 14,012 (78.2%) has N/S ratio of < 1. Similarly, for American goldfinch, we identified 44,732 non-synonymous SNPs and 113,893 synonymous SNPs (N/S ~ 0.39). Out of 16,025 genes, 10,874 (67.8%) has N/S ratio of < 1. We further analyzed the site frequency distribution of SNPs annotated as UTR, synonymous and non-synonymous ( Supplementary Figure 5 ). The distribution showed a standard pattern for a randomly mating natural population, hence we believe this data was suitable for conducting population genetic diversity analysis.

In black-capped chickadees, the inbreeding coefficient (F) averaged 0.11 for the Kent population (0.057 – 0.211) and 0.10 for the Boston population (0.020 – 0.214). These values were normally distributed for Kent and Boston populations and were not significantly different from each other (p = 0.87) ( Supplementary Table 7 ). Similarly, F scores for American goldfinches averaged 0.06 in the Kent population (0.035 – 0.084) and 0.10 (0.045 – 0.167) in the Boston population. These values were also normally distributed and not significantly different from each other (p = 0.13).

The average Weir and Cockerham F_ST estimates across all SNPs loci was ~ 0 between Kent and Boston populations in both species. We further selected SNPs with an F_ST > 0.8 as candidates that represented genomic regions showing significant genetic divergence between geographically distant populations. There were 115 SNPs with an F_ST > 0.8 in black-capped chickadees and 38 SNPs with an F_ST > 0.8 in American goldfinches. Of these SNPs, 89 fell within a 1kb region upstream or downstream of an annotated gene in black-capped chickadees and 21 fell within a 1kb region upstream or downstream of an annotated gene in American goldfinches.

When comparing the list of genes identified as differentially expressed among geographically distinct populations in both species (Kent and Boston) and those possessing highly divergent candidate SNPs (F_ST > 0.8), no genes overlapped. However, we further used these two sets of genes to create a protein-protein interaction network using the STRING database. For black-capped chickadees, there were 174 nodes (proteins representing DEGs and highly divergent genes) and 75 edges (evidence of protein-protein interaction) ( Figure 6 ). For American goldfinches, there were 28 nodes and only 1 edge ( Supplementary Figure 4 ).

The use of wild birds may be necessary when studying species-specific behaviors, ecological interactions, or adaptations that cannot be adequately replicated using alternative models. Wild birds provide valuable insights into natural processes and behaviors that cannot be fully captured in laboratory settings or with alternative experimental approaches. Therefore, replacing wild birds with alternative models may not be feasible without compromising the scientific validity or relevance of the study. In this study, the number of wild birds sampled was minimized using appropriate statistical analysis techniques to ensure that meaningful conclusions can be drawn while minimizing the number of individuals required. Studying gene expression patterns associated with cognitive mechanisms associated with behavior or adaptation necessitated the analysis of brain tissues, and requirement of sacrificing a limited set of birds. This research question could not be answered using alternative methods or non-invasive sampling techniques. The usage of wild birds in this study strictly adhered to the principles of the 3Rs: Replacement, Reduction, and Refinement.

A high contig N50 statistic, over 95% of the original RNA-seq read pairs mapping back to the transcriptome assemblies in both species supported that the reference transcriptomes were constructed successfully. The transcriptome completeness assessment using BUSCO indicated some evidence of missing orthologs; however, this is explained by the fact that our assemblies are derived from only two tissues (brain and liver), and not all genes are expected to be expressed in every single tissue of an organism (Manni et al., 2021). BUSCO analysis also indicated a high degree of duplication, which can be explained due to the nature of the transcriptome that was assembled consisting of multiple isoforms of the same genes.

The number of differentially expressed genes (DEGs) between the Kent and Boston populations was higher in black-capped chickadees compared to American goldfinches. Several factors could be contributing to this result. First, black-capped chickadees and American goldfinches might exhibit different levels of genetic diversity between their geographically separated populations. If black-capped chickadee populations have higher genetic divergence, it would lead to a higher number of DEGs, as this genetic variation may drive gene expression differences. They are also different species with distinct biological and ecological characteristics. The specific genes that drive differences between the two populations in both species may be functionally related to local adaptation, which may be more pronounced in one species than the other. These observations will require further investigation to understand the underlying biological mechanisms driving these differences between the two species.

In black-capped chickadees, among the list of DEGs, we identified functional enrichment for pathways related to nervous system development, such as generation of neurons, neurogenesis, neuron differentiation, neuron projection development and morphogenesis, neuroblast division, and axon development. Overall, the associated genes are known to be involved in both short-term synapse formation, adult brain maintenance, and embryonic and postnatal brain development. Interestingly, some of our DEGs have been identified in previous studies of brain variation in populations of related species of chickadees. Genes upregulated in Kent populations such as PAK3, KALRN, CEP85L, LEF1, FGF13, and GATA3 are known to be associated with cognitive ability in mountain chickadees (Poecile gambeli) (Branch et al., 2022). Similarly, genes upregulated in Boston populations such as NRN1 and CLMN are known to be associated with an improved cognitive phenotype (Branch et al., 2022). Another study by Pravosudov et al. (2013) investigating gene expression in the brain associated with environmental harshness in black-capped chickadees from Kansas and Alaska found a gene we identified to be upregulated in Boston, PLXNA2, to be upregulated in Kansas. However, just because we have identified some of the same genes, does not mean we can definitively say they are associated with environmentally driven variation in cognition or spatial learning and additional measurements of associated trait data in either population will be required in the future.

The DEGs we have identified may be associated with long-term genetic divergence (due to selection-mediated local adaptation or non-adaptive processes such as genetic drift) or plastic responses to short-term environmental differences between two locations. When assessing the environmental variables present in Kent and Boston, overall, they have a similar long-term climate, though they have a great degree of interannual variability. During the month before sampling, the climate in Kent was markedly harsher with greater snow cover and colder temperatures than in Boston. In species of birds that rely on food-caching to survive winter, the harsher an environment is, the greater that bird’s survival depends on accurate retrieval from caches, even short-term caches (Chancellor et al., 2011). The upregulation of genes associated with activity-dependent plasticity and spatial memory in birds that were exposed to colder temperatures and higher snow cover could be a result of this need for enhanced cognition to cope with challenges posed by harsh environmental conditions.

Beyond underlying genetic structure, it is also possible that early developmental conditions rather than environmental fluctuations in these bird’s adult life led to differences in gene expression patterns (Hoshooley et al., 2007). A comparative study examining the epigenetic profiles of these populations could provide valuable insight into the impact of short-term environmental fluctuations on gene expression. Such studies can help us understand the interplay between gene expression, epigenetics, and environmental factors in shaping the cognitive and behavioral traits of these populations. In addition, as food availability is a crucial factor influencing the stress levels of birds during winter, surveying the amount of supplemental feeding or general food availability in the sampling area will be an important consideration for informing and contextualizing future studies on the impact of environmental fluctuations on cognitive and behavioral traits in these birds.

By calculating transcriptome-wide patterns of heterozygosity, we assessed the genetic diversity within and between populations. Two-sample, unpaired t-tests showed that the inbreeding coefficient (F score) between the two black-capped chickadee populations and the two American goldfinch populations in Kent and Boston were not significantly different from each other. Across all individuals, the F value was relatively low. Overall, it seems there is not much inbreeding occurring within these populations as the overall number of homozygous sites is close to the expected number of homozygous sites. This shows there is sufficient standing genetic variation within both populations, indicating a large effective population size in both species. However, a larger sample size would be needed in future studies to get more comprehensive population-wide genetic estimates.

The mean Weir and Cockerham F_ST estimates that compared genetic divergence between populations were also low overall between Kent and Boston populations in both species. This low value of F_ST indicated there is relatively low genetic differentiation in the protein-coding regions between the two populations. A low F_ST value in protein-coding regions indicates sequence conservation and a lack of divergence in functionally important regions (transcribed genes). Genetic conservation of protein-coding regions often suggests that the genes perform vital functions, and any significant changes or mutations in these regions may have deleterious effects on the organism’s fitness (Tatarinova et al., 2016). Studying the whole genome, including non-coding regions under less evolutionary constraint, will provide a more comprehensive and accurate indication of overall genetic divergence between these populations.

Yet, despite a low mean F_ST estimate across all sites, there were still SNPs with high F_ST values indicating the fixation of different alleles between these populations. Between Boston and Kent populations, the American goldfinch brain transcriptome had 762 SNPs with an F_ST above 0.6, 252 with an F_ST above 0.7, 38 with an F_ST above 0.8, 9 with an F_ST above 0.9, and 1 with an F_ST = 1. The black-capped chickadee brain transcriptome had 1,057 SNPs with an F_ST above 0.6, 424 with an F_ST above 0.7, 115 with an F_ST above 0.8, and 14 with an F_ST above 0.9, all of which have an F_ST = 1. Similar to the results from the differential expression analysis conducted, it appears that black-capped chickadee populations from Kent and Boston are genetically more diverged compared to American goldfinches. While genetic differentiation among these geographically separated populations might arise from environmental stressors driving the selection of distinct alleles in each location, it could also result from demographic processes like genetic drift. Additional investigation such as direct assessments of the fitness of these individuals is required to ascertain whether these divergences in coding sequences are under selection or random genetic drift.

Interestingly, none of the SNPs with a high F_ST between populations in both species were found in the vicinity of the genes identified to be differentially expressed in these populations. Because of this, we can reject the hypothesis that sequence variations occurring within the protein-coding regions of the differentially expressed genes are regulating changes in their gene expression levels. However, as RNA-Seq data only spans transcribed regions of the genome, our SNP dataset was not able to capture any mutations occurring in cis-regulatory elements such as promoters, enhancers, and silencers which lie in non-coding regions up and downstream of the target genes. Therefore, we cannot reject or accept the hypothesis that differential expression levels in these candidate genes are due to mutations in the cis-regulatory regions of the genome.

Although our findings revealed no common genes displaying genetic differentiation and alteration in their expression patterns between geographically distant populations, it remains possible that these two sets of genes could still interact through shared biological pathways, or even by functioning as trans-acting factors. By examining the protein-protein interaction networks between DEGs and the genes falling within a 1kb vicinity of the divergent SNPs (F_ST > 0.8) in black-capped chickadees, we identified possible trans-acting factors that may be regulating genes that were differentially expressed. We identified ten transcription factors (BHLHE22, FOSL2, GATA3, ID4, IRX1, IRX6, LEF1, NACA, SOX14, SP5) among the list of differentially expressed genes in black-capped chickadees. Seven were upregulated in Kent (GATA3, ID4, IRX1, IRX6, LEF1, SOX14, SP5) and three were upregulated in Boston (BHLHE22, FOSL2, NACA). These genes may function as trans-regulators of other DEGs. Though they did not have any divergent SNPs within their coding regions, they may still have sequence variation in a non-coding cis-regulatory region, that possibly influences the expression of their downstream target genes.

Three transcription factors (IRX1, IRX6, and LEF1) were found to have protein-protein interactions with one gene containing a divergent SNP, Sonic Hedgehog (SHH). In mammals, SHH plays a large role in both embryonic development and in regulating continual differentiation and proliferation of neuronal precursor cells in adult mammals (Lai et al., 2003; Banerjee et al., 2005). However, the role of SHH in adult neurogenesis in birds has yet to be confirmed. Besides the transcription factors, SHH shares an edge with two other DEGs, FGF13 and OFD1. Though SHH is not a transcription factor, it is a signaling molecule that can regulate other transcription factors (Sasai et al., 2019) and it may be involved in the trans-regulation of the DEGs we have identified.

We have also identified divergent SNPs in the vicinity of three additional transcription factors (TAF1, PURA, and ZF507). TAF1 is a TATA-box binding protein-associated factor, which plays an important role in transcription initiation. Within the brain, TAF1 plays an important role in embryonic neurodevelopmental processes, and mutations have been linked to intellectual disabilities (Grudmundsson et al., 2019). We have identified two DEGs (RERE and POLR2G) having protein-protein interactions with TAF1. Both RERE and POLR2G were upregulated in Kent, hence, sequence variation in the vicinity of TAF1 may lead this gene to have a trans-regulatory effect associated with the upregulation of RERE and POLR2G in Kent populations.

Overall, these results indicate that the genetic architecture underlying the differential gene expression patterns is likely regulatory rather than protein-coding genetic divergence. First, no differentially expressed genes had highly divergent SNPs in their coding regions. Also, some of the DEGs themselves are transcription factors that may play a role in the trans-regulation of other DEGs. Some of the genes within the vicinity of highly divergent SNPs are transcription factors or regulatory molecules that may play a role in the trans-regulation of target DEGs. Additionally, there may be cis-regulator mutations within regulatory regions upstream and downstream of the DEGs themselves. However, though we have identified putative genes that might be involved in the trans-regulation of gene expression, functional studies using gene knockouts would be required to determine the effect they are having in the avian brain. Additionally, in both species, beyond cis- or trans-regulatory control, the differences in gene expression could also be regulated by epigenetic modifications. Further experiments would need to be conducted, such as Whole Genome Bisulfite Sequencing, CHIP-sequencing, or ATAC-sequencing, to examine variation in epigenetic markers such as DNA methylation or histone modification, respectively, between the two populations that may be associated with regulation of gene transcription and subsequent expression levels in either population.