Anterior Pituitary Transcriptome Suggests Differences in ACTH Release in Tame and Aggressive Foxes.

Domesticated species exhibit a suite of behavioral, endocrinological, and morphological changes referred to as "domestication syndrome." These changes may include a reduction in reactivity of the hypothalamic-pituitary-adrenal (HPA) axis and specifically reduced adrenocorticotropic hormone release from the anterior pituitary. To investigate the biological mechanisms targeted during domestication, we investigated gene expression in the pituitaries of experimentally domesticated foxes (Vulpes vulpes). RNA was sequenced from the anterior pituitary of six foxes selectively bred for tameness ("tame foxes") and six foxes selectively bred for aggression ("aggressive foxes"). Expression, splicing, and network differences identified between the two lines indicated the importance of genes related to regulation of exocytosis, specifically mediated by cAMP, organization of pseudopodia, and cell motility. These findings provide new insights into biological mechanisms that may have been targeted when these lines of foxes were selected for behavior and suggest new directions for research into HPA axis regulation and the biological underpinnings of domestication.

Domesticated species exhibit greatly reduced fearfulness, increased social tolerance, and increased resistance to stress compared to their ancestral wild species (Price 2008; Campler ; Driscoll ). These behavioral changes are closely linked to a suite of physiological and morphological changes, including reduced reactivity of the hormonal stress response, white spotting, curling of the tail, and shortening of the cranium, which have together been referred to as the “domestication syndrome” (Hare ; Wilkins ). The closely correlated behavioral, endocrinological, and morphological changes across multiple species suggest that perturbations at the level of networks of genes, working together for specific biological effects, may be implicated in the domestication process. Therefore, we expect that understanding the biological basis of the endocrinological changes in one domesticated species may shed light on the mechanisms targeted during the domestication process more generally, and potentially also on the biological underpinnings of variation in resilience to stress in animals and humans (Hare ).

The characteristic reduction in fearfulness of domesticated animals has been closely linked to reduced reactivity of the hypothalamic-pituitary-adrenal (HPA) axis—the hormonal cascade associated with the stress response in mammals (Trut ; Hare ; Wilkins ). Specifically, reductions in circulating levels of adrenocorticotrophic hormone (ACTH), released by the anterior pituitary, have been demonstrated in several domesticated species (Gulevich ; Trut ; Kaiser ). ACTH stimulates release of glucocorticoids by the adrenals, and reduction in glucocorticoids in domesticated species has also been shown (Albert ; Trut ; Kaiser ). These differences in domesticated and wild HPA axis function have been explored in multiple animal models. Domesticated guinea pigs, in response to psychosocial stressors, have adrenal glucocorticoid responses both decreased in magnitude and returning to baseline after a shorter time (Kaiser ). Similar findings have been demonstrated with experimentally domesticated Norway rats, which, after many generations of selection for tameness, display markedly reduced aggressiveness to humans as compared to wild rats as well as reduced glucocorticoid levels (Naumenko ; Albert ; Prasolova ) and ACTH levels (Shikhevich ). However, the mechanisms driving the difference in HPA activity between these tame animals and their wild cousins are unknown.

In perhaps the best-known domestication experiment, farmed foxes (Vulpes vulpes) have been selectively bred for >50 generations to recapitulate the dog domestication process. One population of foxes was selected for lack of fear of, and for sociability toward, humans, while a second population of foxes was selected for aggression toward humans. A comparison of these tame and aggressive fox populations demonstrates similar changes in the HPA axis as those shown between domesticated guinea pigs and their wild ancestors, as well as between tame and aggressive rats (Trut , 2009). Compared to the aggressive line, the tame foxes display reduced basal ACTH levels (Gulevich ; Trut ), decreased ACTH response to stress, and shortened time to return to baseline ACTH levels (Oskina and Plyusnina 2000), as well as markedly reduced glucocorticoid elevation and shortened time to return to baseline glucocorticoid levels after a stressor (Trut ).

Although the hypothalamus is classically considered the “top” of the HPA axis, the pituitary, situated to amplify signals from the hypothalamus before release to the systemic circulation, is sometimes known as the “master gland” (Melmed 2010) and may contribute significantly to regulation of the HPA axis. Inhibition of the HPA hormonal cascade at the level of the pituitary may occur in at least three ways, any or all of which could have been selected for during the domestication process: reduced synthesis of ACTH, reduced secretion of ACTH, and reduction of the number of cells that synthesize ACTH (Rizzoti 2010). An examination of the details of each of these processes will provide predictions for how they might vary in wild compared to domesticated pituitaries (i.e., the pituitaries of wild vs. domesticated animals). ACTH synthesis and secretion are triggered by the binding of corticotrophin releasing hormone (CRH) to CRH receptors in the anterior pituitary (Grammatopoulos 2012). ACTH synthesis begins with the production of the proopiomelanocortin (POMC) precursor molecule, which is regulated by a suite of transcription factors (Drouin 2016) and is cleaved into ACTH by proprotein convertases (Benjannet ), which may also be thought to regulate ACTH production. Therefore, the domesticated pituitary might contain different concentrations of CRH receptors, POMC, POMC-stimulating transcription factors, or POMC convertases compared to the aggressive pituitary.

Multiple mechanisms regulate ACTH secretion through exocytosis, but the main signal for release of hormone-filled secretory vesicles is a rise in calcium from both intra- and extracellular sources, and this rise in calcium is primarily regulated by cAMP (Burgoyne and Morgan 2003). Some pituitary hormone-releasing cells may also facilitate the release of hormones into the systemic circulation by extending pseudopodia and becoming motile to contact blood vessels (Navratil ). Additionally, pituitary endocrine cells are known to form networks, potentially via gap junctions, to coordinate hormone release (Hodson ,b). Therefore, the domesticated pituitary may demonstrate different levels of motility, cellular adhesion and communication, or cAMP activity as compared to the wild pituitary.

Finally, tame and aggressive pituitaries may display different population sizes of ACTH-secreting cells (corticotrophs), with a larger population associated with greater ACTH production and/or release. Stem cells in the pituitary differentiate into specialized hormone-producing cells, including corticotrophs, throughout life at a low rate (Florio 2011). Therefore, key transcription factors involved in corticotroph cell fate determination may differ in the wild and domesticated pituitary.

To identify mechanisms involved in the regulation of reduced reactivity of the HPA in domesticated animals, we analyzed gene expression in anterior pituitary of tame and aggressive farmed foxes. Previously, tame fox pituitaries have been shown to have comparable amounts of POMC mRNA and ACTH protein to aggressive fox pituitaries (Gulevich ), despite decreased circulating concentrations of ACTH. Therefore, we specifically investigated whether the changes in ACTH release in tame foxes are associated with ACTH secretion (exocytosis) and/or corticotroph differentiation, rather than ACTH synthesis from POMC mRNA. Characterization of genes and gene networks with different anterior pituitary activity in tame and aggressive foxes will provide candidate mechanisms for differences in HPA axis regulation between domesticated and wild animals and may inform future studies into the biological underpinnings of domesticated animal behavior, such as reduced fearfulness, increased social tolerance, and increased resistance to stress.

Materials and Methods

Animals and sample preparation

Animals were maintained at the experimental farm of the Institute of Cytology and Genetics (ICG) in Novosibirsk, Russia. All animal procedures at the ICG complied with standards for humane care and use of laboratory animals by foreign institutions. The study was approved by the Institutional Animal Care and Use Committees (IACUC) of the University of Illinois at Urbana-Champaign. Anterior pituitary tissue was dissected from six tame and six aggressive sexually naive 1.5-yr-old males. The foxes were selected based on their relatedness and behavioral scores. Most of the selected foxes did not share parents and grandparents, with the exception of two aggressive foxes that had one common relative (a parent of one fox and a grandparent of another fox), and two pairs of tame foxes, each of which shared one grandparent. All of the selected foxes from the tame population belonged to the “elite” group of domesticated foxes, i.e., had the highest behavioral scores, while all selected foxes from aggressive population had the lowest (most aggressive) behavioral scores (Trut ; Kukekova ). Foxes were killed using sodium thiopental; immediately afterward, the brain case was opened with a saw and the brain was removed, leaving the pituitaries in situ. The anterior and posterior pituitaries were immediately dissected out and placed into RNAlater (Life Technologies, Grand Island, NY). The samples were stored at −70°. Total mRNA was extracted using RNeasy Mini Kit (Qiagen, Valencia, CA), according to the manufacturer’s protocol.

Sequencing and quality analysis

One microgram of high quality RNA from each sample was used for sequencing. Stranded RNA-seq libraries were prepared using TruSeq SBS Sequencing kit version 3 (Illumina, San Diego, CA). Libraries were barcoded and pooled and sequenced on two lanes on a HiSeq2500. Reads were single-end, stranded, and 100 nt in length. Sequencing results were processed by CASAVA 1.8 (Illumina, San Diego, CA).

Data quality, including base quality per position across reads, GC content, and distribution of sequence length, was initially assessed with FastQC (Andrews 2010). Reads were processed with flexbar (Dodt ) in two passes: the first to trim adapters, remove low quality reads, and remove reads <35 bp in length, and the second to remove polyA tails. Subsequently, reads that mapped to fox mitochondrial DNA sequences from NCBI (accession numbers JN711443.1, GQ374180.1, NC_008434.1, and AM181037.1) using Bowtie2 (Langmead and Salzberg 2012) were discarded. Similarly, any remaining reads that mapped to ribosomal DNA sequences were discarded.

Evaluation and validation of differential expression

Reads were aligned to Canis familiaris 3.1 (CanFam3.1) with Ensembl 1.79 annotation using TopHat 2.0.13 (Kim ). The Ensembl annotation used is available at . The fox and the dog ancestor, the gray wolf, diverged ∼7–10 million yr ago (Wayne 1993). Therefore, to account for expected variation between our fox sequences and the CanFam3.1 reference genome, we added the following TopHat2 parameters to the command line, to allow a more robust cross-species alignment with the dog genome: -N 3 –read-edit-dist = 3. These parameters increased the number of nucleotide mismatches and maximum edit distance allowed per read. Gene expression levels were assessed using HTSeq-count 0.6.1p1 (Anders ). Differentially expressed genes between samples from tame vs. aggressive foxes were evaluated using DESeq2 1.10.1 (Love ), and results were corrected for multiple testing using the Benjamini-Hotchberg algorithm. Normalized gene counts (FPKMs) were assessed with CuffNorm 2.2.1 (Trapnell ) using quartile normalization.

RNA from all 12 samples used for sequencing was also used for quantitative real-time PCR (RT-qPCR). Gene assays used were AKAP7, ANKRD6, BRCA2, CDH4, ITGA1, LAMA2, NELL2, PDE7B, and SPOCK1. Primers for qPCR (Table S1 in Supplemental Material, File S1), placed in adjacent exons and spanning exon-intron boundaries where possible, were designed with Primer3 (Koressaar and Remm 2007; Untergasser ). cDNA was synthesized from 1 μg of RNA with random hexamers using the ThermoScript RT-PCR Kit (Invitrogen, Carlsbad, CA). RT-qPCR was performed on 96-well plates using an Applied Biosystems StepOnePlus Real Time PCR System (Applied Biosystems, Foster City, CA). Wells were loaded with 10 µl solution containing 5 µl SYBR Green reagent (Applied Biosystems), 0.1 µg cDNA, and either 2 mol forward and 2 mol reverse primers for the control gene or 0.6 mol forward and 0.6 mol reverse primers for the test genes. All samples were analyzed in triplicate. Initial denaturation was performed at 95° for 20 min, followed by 40 cycles of 95° for 3 sec and 60° for 3 sec, followed by a melting curve analysis at 60–95°. Gene expression was normalized to the expression of SDHA—a gene that did not show differential expression between tame and aggressive foxes and that had comparable expression between all foxes in the RNA sequencing results. Mean ΔΔ-Ct values for tame foxes were compared to mean ΔΔ-Ct values for aggressive foxes for each gene using the Mann-Whitney test.

Functional analysis

Functional analysis was performed for each tissue with ClueGO version 2.2.5 (Bindea ); Gene Ontology (GO) files were downloaded on May 25, 2015. The ClueGO evaluation was performed against a background list of genes with a mean expression level in pooled tame and aggressive pituitary samples reported to be >0 by DESeq2. Resulting p-values were adjusted using the Benjamini-Hochberg correction. Enriched GO terms were reviewed for association with stimulus of hormone release, POMC regulation and processing, and corticotroph differentiation.

Weighted gene coexpression network analysis

Modules of genes with highly correlated expression patterns were sought using weighted gene coexpression network analysis. We expect these modules to correspond to networks of genes that interact in shared biological processes (Langfelder and Horvath 2008). The full set of genes expressed in the pituitary reads was filtered to remove genes with <10 FPKMs in each of three or more samples. We then constructed unsigned weighted gene coexpression modules using the WGCNA package in R (Langfelder and Horvath 2008). The blockwiseModules function was run with the Pearson correlation coefficient and a soft thresholding power of 12. The resulting gene groups, or modules, were named by assigning them arbitrary colors. The first principle component of each module was designated its eigengene. Eigengenes for each module were divided into samples from tame and aggressive foxes; these two groups were compared with a two-tailed Student’s t-test and their p-values were adjusted using the Benjamini-Hochberg correction. Modules selected for investigation were analyzed for gene enrichment using ClueGO as described above. Additionally, highly connected “hub” genes in these modules were identified using Cytoscape (Shannon ) to select the genes with the largest number of connections to other genes in the same network, with the expectation that these “hub” genes would represent regulatory genes with control over many genes in that module.

Alternative splicing analysis

Differences in alternative splicing frequencies between tame and aggressive pituitary reads were explored with rMATS (Shen ). Because rMATS requires reads of equal lengths, postfiltered reads were further filtered using SAMtools 1.1 (Li ) to remove all reads <100 nucleotides in length. BAM files (previously aligned to CanFam3.1 for differential gene expression calling) and Ensembl annotation, version 1.79, were used with rMATS to investigate reads crossing splicing junctions. Differences in frequencies of reads including, or skipping, specific exons were used to identify exons that were skipped significantly more frequently in tame or aggressive pituitary reads. P-values were adjusted using Benjamini-Hotchberg. The list of differentially skipped exons was filtered to remove exons covered by <10 reads in either tame or aggressive samples. Differences in skipping were reported as mean percentage skipped exons, in other words, (mean number of reads skipping the exon)/(mean number of reads skipping the exon + mean number of reads including the exon). Exons identified as significantly more frequently skipped in tame or aggressive samples were further investigated by searching for their amino acid sequence using their nucleotide sequence in blastx () to identify isoforms that included, or did not include, the exon in question. Conserved domains that the exon in question fell within were identified using NCBI’s Conserved Domains Database (https://www.ncbi.nlm.nih.gov/cdd/).

Data availability

The primary sequencing data have been deposited in the GenBank Sequence Read Archive (SRA) under accession number SRP107226.

Results

Gene expression in the anterior pituitary of tame and aggressive foxes

Sequencing of fox pituitaries produced 206 million reads total for six tame samples (194 million reads after filtering) and 200 million reads total for six aggressive samples (191 million reads after filtering) (Table S2 in File S1). The filtered fox reads were mapped to the Canis familiaris 3.1 reference with a mean alignment rate of 78.98% (Table S3 in File S1). The rate of multiple alignments was low, with a mean of 2.64%. Analysis of filtered pituitary reads identified 18,940 distinct genes expressed in this tissue (Table S4). Analysis of gene expression level identified 11 genes with notably higher expression (Figure 1). A GO analysis was performed on the 0.1% most highly expressed genes (180 genes). Of the most significantly enriched terms [false discovery rate (FDR) < 10.0e−9; nine terms], eight were related to RNA processing or splicing (Table S5). One of these genes, EEF1A1, was significantly more highly expressed in the tame than the aggressive fox pituitary. The most highly expressed gene, by an order of magnitude, was POMC, which codes for the precursor protein for ACTH and other peptides released by the pituitary.

Figure 1

Most highly (0.01%) expressed genes in fox pituitary by normalized gene count, assessed by DESeq2. The most highly expressed gene, POMC, is excluded because its order of magnitude larger expression level than the next largest gene results in a scale that does not allow detailed values to be displayed for genes with smaller expression levels.

DESeq2 analysis of filtered pituitary reads identified 346 differentially expressed genes (DE genes) (FDR < 0.05); of these, 191 (55%) were upregulated in aggressive pituitary, and 155 (45%) were upregulated in tame pituitary. As predicted by Gulevich ), POMC was not differentially expressed (expression was slightly greater in the tame samples, with log2 fold change = 0.003, FDR = 0.996). Log2 fold changes among DE genes ranged from −0.81 to 0.84 (mean −0.043, SD 0.37) (Table S6). Clustering analysis of DE genes clearly differentiated tame and aggressive samples (Figure S1 in File S1). Additional visualization using a volcano plot showed that the distribution of log2 fold changes and significance was as expected, and principal components (PC) analysis grouped tame samples separately from aggressive samples (Figure 2 and Figure S2 in File S1).

Figure 2

Volcano plot of gene expression in fox pituitary. Red dots represent genes significantly differentially expressed in tame vs. aggressive fox lines at FDR < 0.05. Green dots represent genes that are both significantly differentially expressed and have a log2 fold change >0.5. Outlier genes are labeled with their gene symbols.

Nine DE genes (AKAP7, ANKRD6, BRCA2, CDH4, ITGA1, LAMA2, NELL2, PDE7B, and SPOCK1), with a range of FDR values and log2 fold changes, were selected for validation with quantitative PCR (qPCR). Analysis with RNA-seq data indicated that four of these genes were upregulated in aggressive foxes and five in tame foxes. Comparison of RNA-seq and qPCR results found one gene, ITGA1, which was shown to be upregulated in tame foxes according to the RNA-seq results but downregulated in tame foxes according to the qPCR results, though this result was not statistically significant. However, this gene had a very small log2 fold change in the DESeq2 output (0.21). The remaining eight of the tested genes showed expression differences in the same direction in both RNA-seq and qPCR results (see Figures S3 and S4 and Table S7 in File S1). The smallest P value (P = 0.06) in comparison of qPCR values for tame and aggressive samples was obtained for SPOCK1—the gene with the largest log2 fold change and the second smallest FDR (according to RNA-seq results) tested.

ClueGO analysis:

ClueGO analysis was performed on the 328 DEG that had an associated gene symbol; a background list was used of 18,940 genes for which DESeq2 demonstrated expression in anterior pituitary. The DEG were associated with 394 GO terms, of which 74 GO terms were significantly enriched (FDR < 0.05) (Table S8).

The most significantly enriched terms in the GO analysis at FDR < 0.01 (Table 1) cluster into five groups with 11 terms total. Three of these groups are represented by a single GO term each (pseudopodium organization, regulation of dendritic cell differentiation, and granulocyte chemotaxis). The most significantly enriched GO term relates to pseudopodia—a generic term for different types of membrane projections known to function in cellular motility, as well as in release of lutenizing hormone from anterior pituitary gonadotrophs (Navratil ). This GO term includes the gene CDC42EP2, which induces pseudopodia formation (Hirsch ); SRGAP2 and SRGAP3, which regulate cell migration through modification of cell protrusion morphology (Guerrier ; Endris ); and WASH1, which promotes pseudopod extension (Zech ).

Table 1

GO terms enriched in genes differentially expressed in tame vs. aggressive pituitary tissues at FDR < 0.01

Group	GO ID	Term Description	Adjusted P-Value	DEG Annotated with this Term	% Genes in this Term Differentially Expressed
1	GO:0031268	Pseudopodium organization	0.007	CDC42EP2, SRGAP2, SRGAP3, WASH1	28.6
2	GO:0007259	JAK-STAT cascade	0.007	ASPN, CISH, ERBB4, FLRT1, FLRT3, IL6ST, LRRC4B, PTGER4, SOCS2, STAT3a, TNFRFSF1A	8.3
	GO:0097696	STAT cascade	0.007		8.3
	GO:0046425	Regulation of JAK-STAT cascade	0.009		9.2
	GO:1904892	Regulation of STAT cascade	0.009		9.2
3	GO:2001198	Regulation of dendritic cell differentiation	0.008	LGALS3, LGALS9, TMEM176b	50
4	GO:0006614	SRP-dependent cotranslational protein targeting to membrane	0.008	RPL17, RPL23A, RPL26, RPL4, RPL7A, RPS24, RPS3A, ZFAND2Bb	10.7
	GO:0006613	Cotranslational protein targeting to membrane	0.008		10.5
	GO:0045047	Protein targeting to ER	0.008		10.4
	GO:0016259	Selenocysteine metabolic process	0.009		11.7
5	GO:0071621	Granulocyte chemotaxis	0.010	IL17RA, IL17RC, ITGA1, LGALS3, NCKAP1L, PDE4B, TGFB2, THBS4, VAV3	10.1

STAT3 is not associated with GO:1904892.

ZFAND2B is not associated with GO:0016259.

The second of the three single-term GO groups at FDR < 0.01 is regulation of dendritic cell differentiation, which includes genes for two galectins related to cellular motility [LGALS3 and LGALS9 (Sano ; Bi )] and a gene related to the cross-presentation of antigen [TMEM176B (Segovia )]. The last term of the single-term GO groups, granulocyte chemotaxis, contains two genes for receptors for interleukin 17A (IL-17), a pro-inflammatory cytokine, IL17RA and IL17RC (Gaffen 2009). This term also contains a galectin, LGALS3; an integrin, ITGA1; VAV3, a gene associated with cell migration and integrin signaling (Sachdev ; Sindrilaru ) and implicated in the formation of lamellipodia and membrane ruffles (Movilla and Bustelo 1999); PDE4B, a phosphodiesterase, related to cAMP regulation (Houslay and Adams 2003); and THBS4, related to cellular adhesion, migration, and attachment (Adams 2001; Mosher and Adams 2012).

Four terms among the most significant GO terms are related to the JAK-STAT pathway (JAK-STAT cascade, STAT cascade, regulation of JAK-STAT cascade, and regulation of STAT cascade). The JAK-STAT pathway is associated with the interleukin immune response (Subramaniam ; Sansone and Bromberg 2012), cell migration (Pocha and Montell 2014), and the cAMP pathway (Chesnokova and Melmed 2002). Four other related GO terms from Table 1 (SRP-dependent cotranslational protein targeting to membrane, cotranslational protein targeting to membrane, protein targeting, and selenocysteine metabolic process) are related to protein targeting and signal recognition proteins (SRP). Protein targeting is a process by which proteins are translocated to an organelle identified by SRPs during translation.

GO analysis was also performed separately on genes significantly upregulated in the tame samples, and genes significantly upregulated in the aggressive samples. Three GO terms were significantly enriched for genes upregulated in the aggressive samples (Table S9); these terms had an overlapping set of DEG associated with them, including many of the genes in group 1 in Table 1. These three terms all related to pseudopodia. Twenty-seven GO terms were significantly enriched among the genes upregulated in tame samples (Table S10). Many of these terms had overlapping DEG, and could be sorted into four groups, containing genes associated with: white blood cell migration (containing genes overlapping with groups 3 and 5 in Table 1); protein targeting to the endoplasmic reticulum (containing genes overlapping with group 4 in Table 1); transcription and translation (also containing genes overlapping with group 4 in Table 1); and the humoral immune and complement response. This final group of terms contained genes not identified in any of the groups in Table 1, except for LGALS9, from group 3. The other genes in this group were BCL3, C1R, C1S, C7, CFI, PAXIP3, and SLC11A1.

Exploration of additional terms that were significantly enriched in DEG (Table S8) identified terms associated with cell migration (including pseudopodia and cytoskeletal organization), cell adhesion, exocytosis (including formation of exocytotic vesicles), cAMP regulation, and calcium regulation. Specific terms include: neutrophil migration, negative regulation of neuron migration, granulocyte migration, pseudopodium assembly, cell motility involved in cerebral cortex radial glia guided migration, neutrophil chemotaxis, myeloid leukocyte migration, regulation of neuron migration, leukocyte chemotaxis, and glial cell migration. Several terms are also included associated with cellular adhesion and cellular interactions with the extracellular matrix, critical components of cellular migration: extracellular structure organization, regulation of integrin activation, extracellular matrix organization, integrin activation, and positive regulation of cell adhesion mediated by integrin. Finally, several terms relate to biosynthesis of cAMP—an important signaling molecule in ACTH release by exocytosis (Reisine ): cAMP metabolic process, cAMP biosynthetic process, positive regulation of cAMP biosynthetic process, positive regulation of cAMP metabolic process, and regulation of cAMP biosynthetic process.

As part of an investigation of the hypothesis that changes in the tame fox pituitary would primarily relate to hormone release, terms associated with CRH-dependent regulation of POMC transcription, processing of the POMC protein, and differentiation of corticotrophic cells were sought in the list of GO terms enriched for DEG. No terms related to these topics were found to be enriched at FDR < 0.05.

Analysis of genes involved in exocytosis, POMC regulation, and cell differentiation:

DEG related to exocytosis, regulation of POMC expression, and cell differentiation were sought, to investigate which of these processes were more extensively represented in the differences in gene expression between tame and aggressive fox pituitaries (Table 2). The gene with the largest log2 fold change was GRK7 (upregulated in aggressive pituitary), a member of the G-protein coupled receptor kinase family. Genes from this family are known to serve as important regulators of phosphorylation of the CRH receptor, activation of which triggers ACTH synthesis and release (Grammatopoulos 2012). AVP, which codes for arginine vasopressin—a neuropeptide with paracrine activity in the anterior pituitary known to potentiate ACTH release (Liu ) but not synthesis (Levin )—was upregulated in the tame pituitary. DEG associated with membrane proteins of exocytotic vesicles included VAMP1 (Calakos ), and two synaptotagmins [SYT10 and SYT13 (Brose )]. Therefore, tame and aggressive pituitaries demonstrated potential differences in regulation of ACTH release via CRH receptor regulation, vasopressin signaling, and exocytotic vesicle formation, as well as potential differences in ACTH synthesis via CRH receptor regulation.

Table 2

DEG in tame vs. aggressive fox pituitary tissue involved in exocytosis, POMC regulation, and cell differentiation

Gene Name	FDR	Log2 Fold Change
GRK7	1.17e−011	−0.81
NPR3	8.46e−005	−0.59
CALCRL	1.02e−005	−0.45
SYT10	0.04	−0.39
PDE7B	0.02	−0.25
AKAP7	0.002	−0.23
RAPGEF4	0.03	−0.21
VAMP1	0.04	0.28
PDE4B	0.02	0.28
SYT13	0.01	0.33
ADCY7	0.04	0.38
AVP	0.03	0.41
ADCY8	0.01	0.42

Results are ordered by log2 fold change. Negative log2 fold change values indicate genes upregulated in aggressive samples; positive log2 fold change values indicate genes upregulated in tame samples.

Regulation of cAMP was specifically investigated, as cAMP serves as a major cellular signal for exocytosis. DEG associated with cAMP regulation included two members of the adenylyl cyclase family of genes that synthesize cAMP (ADCY7 and ADCY8); two members of the phosphodiesterase family of genes that break down cAMP (PDE4B and PDE7B); G-protein receptors that activate adenylyl cyclase [CALCRL/CRLR (Smith ) and P2RY11 (Qi )]; a G-protein receptor which inhibits adenylyl cyclase [NPR3 (Anand-Srivastava 2005)]; a guanine nucleotide exchange factor that functions in cAMP-dependent exocytosis [RAPGEF4/EPAC2 (Sugawara )]; and a scaffolding protein involved in cAMP compartmentalization [AKAP7/AKAP15 (Kritzer )]. Additionally, the previously discussed DEG AVP is known to act through the cAMP/PKA pathway in the anterior pituitary (Liu ) in its effects on POMC regulation. Therefore, differences in cAMP regulation are represented in multiple ways in this data set, from synthesis to metabolism to compartmentalization.

Notably, genes important in CRH-dependent regulation of POMC transcription, such as NGFIB, NURR1, and NOR1 (Bonfiglio ), were not differentially expressed, nor were genes associated with processing of POMC into ACTH, such as PC1 and PC2 (Benjannet ), CPE (Cool and Loh 1998), and PAM (Ciccotosto ).

Genes associated with pituitary, and, specifically, corticotroph (ACTH-producing), cell differentiation were investigated to determine the likelihood that tame and aggressive pituitaries have different population sizes of corticotrophic cells. No genes directly associated with corticotroph differentiation were identified as significantly differentially expressed. However, TBX19/TPIT, a gene known to be essential in the establishment of the POMC-expressing lineage of corticotrophs and melanotrophs (Pulichino ), was significantly upregulated in aggressive pituitary before adjustment at P = 0.047 (FDR = 0.32). PAX7, a gene known to suppress corticotroph differentiation in favor of melanotroph differentiation, (Budry ), was also significantly upregulated in aggressive pituitary only before adjustment (P = 0.03, FDR = 0.28).

Weighted gene coexpression network analysis results

Weighted gene coexpression network analysis constructed 66 gene modules with activity in fox anterior pituitary. Two modules, Antique White 4 and Pale Turquoise, had significantly (FDR < 0.10) different eigengenes between tame and aggressive samples (Table 3). Both modules had a negative t-statistic, indicating that they are more active in aggressive than tame samples.

Table 3

Modules of coexpressed gene networks with significantly different activation in tame vs. aggressive fox pituitary tissue

Module Name	t-Statistic	FDR	Number of Genes in Module	Percent Genes in Module Differentially Expressed
Antique White 4	−4.2	0.09	43	18.6
Pale Turquoise	−4.1	0.09	80	20.0

The Antique White 4 module had eight DEG in total, including SYT10, a calcium sensor involved in exocytosis of secretory vesicles (Cao ); RNPEP, an aminopeptidase (Piesse ); and two genes associated with apoptosis, CPPED1/CSTP1 (Zhuo ) and TRAF4 (Sax and El-Deiry 2003). (For DEG, see Table S6; for genes in WGCNA modules, see Table S11 and Table S12) The Antique White 4 module had no enriched GO terms, possibly due to the small number of genes (43) it contains. Hub (i.e., highly connected) genes in this module included an aminopeptidase, RNPEP (Piesse ), which was both the most highly connected gene and upregulated in the aggressive pituitary; ZIC5, associated with neural crest cell differentiation (Nakata ; Inoue ); SUPT6H/SPT6, a transcription elongation factor (Diebold ); KLHL2, involved with cell membrane process formation (Williams ); and ATG5, involved with autophagic vesicle formation (Walczak and Martens 2013). Additional genes in this module associated with exocytosis, cell migration, and cell–cell adhesion include ANXA7 (Caohuy ), LYN (Malik ; Nakata ), and RAP1B (Mandell ); both LYN and RAP1B are also associated with integrin signaling (Malik ; Mandell ). Overall, this module demonstrated relationships to processes associated with exocytosis, apoptosis, vesicle formation, cell migration, and cell–cell adhesion.

The Pale Turquoise module had 16 DEG. These included a number of genes associated with exocytosis and migration, such as genes associated with the anchoring or biosynthesis of cAMP [AKAP7 (Kritzer ) and FLNB (del Valle-Pérez )]; AVP, a signaling hormone for ACTH release (Denef 2008); an integrin and a laminin, two proteins known to work together in the formation of pseudopodia [ITGA6 (Mercurio ) and LAMA2 (Gu )]; and a cadherin-associated cell–cell adhesion gene, PTPRT (Besco ). There were two significantly (P < 0.01) enriched GO terms for this module, regulation of cytokinesis and intrinsic apoptotic signaling pathway (Table S13). Hub genes in this module included genes associated with puberty [WDR11 (Kim and Layman 2011)] and carnosine peptidase activity [CNDP1/CN1 (Teufel )]. Hub genes that were also differentially expressed included LAMA2, PTPRT, and WLS [Wnt protein regulation and secretion (Bänziger )]. Therefore, this module also had strong associations to genes involved in exocytosis and migration, as well as apoptosis.

Notably, the Brown module, which showed increased activation in tame samples only before correction (P = 0.03, FDR = 0.20) suggested functionality related to alternative splicing of mRNA. This module was associated with three GO terms significantly enriched for module genes that were related to this function, mRNA 3′-splice site recognition, regulation of RNA splicing, and alternative mRNA splicing, via spliceosome. The Brown module additionally contained multiple hub genes associated with mRNA alternative splicing and regulation, including TSEN2/SEN2, TNKS1BP1/TAB182, and PRCC (Trotta ; Skalsky ; Lau ). The mRNA splicing functionality in this module may be associated with cell motility, as significantly enriched GO terms in the Brown module also included the terms lamellipodium organization, ruffle organization, regulation of fibroblast migration, regulation of focal adhesion assembly, wound healing/spreading of cells, and negative regulation of cell-substrate adhesion. Lamellipodia are pseudopodia formed on the leading edge of migrating cells (Krause and Gautreau 2014); membrane ruffling allows the cell to form a motile surface to network with other cells (Ridley 1994). The Brown module had a positive t-statistic of 2.7, indicating possible greater activation in tame than aggressive samples.

Overall, gene network analysis with WGCNA demonstrated two modules, Antique White 4 and Pale Turquoise, with significantly increased activity in aggressive samples. These modules were associated with exocytosis, cell migration, and cell adhesion; both modules contained genes suggestive of the importance of integrins in these functionalities. Both modules also indicated the importance of apoptosis, suggesting that increased cellular turnover may have some significance in differing pituitary function between the lines. Finally, analysis of the Brown module suggested the importance of alternative splicing of mRNA in anterior pituitary, which may or may not show greater activity in the tame population. For this reason, an alternative splicing analysis was undertaken.

rMATS was used to identify 36 genes with exons skipped at different frequencies in tame compared to aggressive fox pituitary reads at FDR < 0.05 (Table S14), and 12 genes at FDR < 0.001 (Table 4). Exon counts for these 36 genes with significantly differentially skipped exons ranged from 4 to 46 exons total (median 13, SD 9.3). Of the genes with differentially skipped exons, only one [MTHFSD, an RNA-binding protein associated with amyotrophic lateral sclerosis (MacNair )] is itself differentially expressed in pituitary reads. Investigation of GO terms revealed that nine terms enriched for DEG also contained genes with differentially skipped exons. Of these, four terms contained multiple genes with differentially skipped exons: establishment of protein localization to organelle (ABLIM3, SH2D4A, POT1; this term is functionally related to the highly significantly enriched group of four terms related to protein targeting); gliogenesis (ADAM22, NFIB); and extracellular matrix organization and extracellular structure organization (both containing KLKB1 and MFAP5). Nine genes with differentially skipped exons were also included in WGCNA modules associated with the state of being tame, including three genes in the Brown module (DAGLB, MYCBPAP, and STAUI). GO enrichment analysis of the genes with differentially skipped exons revealed no enriched terms.

Table 4

Genes containing exons skipped at different frequencies in tame compared to aggressive fox pituitary reads (FDR < 0.001), including conserved domains overlapping the skipped exon

Gene	Exon Location	Percent Reads Skipping Exon (Tame)	Percent Reads Skipping Exon (Aggr)	FDR	Skipped in Isoforms	In Conserved Domain(s)
MFAP5	chr27:37,034,668-37,034,704	46	0	0	X2, X3	MAGP
SPICE1	chr33:17787918-17788056	13	0	2.75E−10	None	SPICE
MECR	chr2:71398391-71398452	0	9.8	1.39E−07	None	ETR, Qor
SPATA20	chr9:26,472,931-26,473,143	12	1.9	5.57E−06	X4	YyaL
KLKB1	chr16:44493733-44493873	20	0	7.98E−00	None	Tryp_SPc
THYN1	chr5:996686-996755	12	2.3	1.81E−05	None	Eve
FBRSL1	chr26:542234-542303	13	2.8	1.29E−04	None	None
TRIM65	chr9:4716736-4716911	14	2.8	1.40E−04	None	None
TMEM198	chr37:26056452-26057028	2.6	0.31	3.22E−04	None	DUF4203
GANC	chr30:9295384-9295464	10	1.8	3.60E−04	None	GH31_GANC_GANAB_alpha, Glyco_hydro_31
SORBS1	chr28:8892641-8892725	67	93	3.76E−04	X7, X8, X12	None
TCFL5	chr24:46680735-46680982	7.1	1.3	8.40E−04	None	None

Of the genes listed in Table 4, MFAP5/MAGP-2 promotes cell adhesion and spreading via integrin binding (Gibson ), and SORBS1/CAP is involved in integrin-mediated cell adhesion (Ribon ). At FDR < 0.05, several genes related to cell migration, cell–cell adhesion, and calcium signaling contained differentially spliced exons. These included: ADAM22, a disintegrin (D’Abaco ); CFAP97/HMW, associated with cilia and flagella (Soulavie ); EFMA5, integrin-dependent cell-cell adhesion (Davy and Robbins 2000); ENKUR, calcium signaling and flagella (Sutton ); ERG, adhesion and cell morphology change (McLaughlin ); UNC5D, cell-cell adhesion/repulsion and motility (Jarjour ). Exons alternatively spliced at FDR < 0.05 in genes related to RNA processing included ELAVL4 (Ince-Dunn ) and FBRSL1 (Baltz ). Additional interesting genes just above significance level were RPH3AL/NOC2 (P = 0.10), associated with exocytosis (Matsunaga ), and POLK (P = 0.07), a DNA repair polymerase (Ogi and Lehmann 2006).

Discussion

RNA sequencing of pituitary tissue from tame and aggressive foxes identified expression differences in genes related to cAMP regulation, pseudopodia formation, and cell–cell adhesion. These findings suggest that the pituitaries of experimentally domesticated foxes may differ from the pituitaries of a lineage of foxes selected for increased aggression in regulation of exocytosis. Specifically, cAMP signaling for vesicle release may be modified in the tame foxes to reduce the release of ACTH during the stress response. Additionally, tame fox pituitary cells may differ in their motility or their ability to form networks for cell–cell communication, potentially to coordinate ACTH release, as indicated by differences in the expression of genes related to the organization of pseudopodial membrane extensions and cell adhesion. Notably, significant differences in lineage specification markers were not found between the two populations, suggesting that differentiation of pituitary cells involved in the HPA axis did not contribute to the differences in function between the tame and aggressive pituitary. However, nonsignificant differences were seen, and, due to the small relative size of the population of actively differentiating cells in the anterior pituitary (Florio 2011), significant differences may be difficult to detect in a mixed sample of cell populations such as this one. Additionally, one gene coexpression module suggested differences in alternative splicing functionality between tame and aggressive fox lines, and differentially skipped exons were identified in genes associated with cell migration, cell–cell adhesion, and calcium signaling. Therefore, differential splicing in the tame fox pituitary may affect similar functionality as do gene expression differences.

Alignment and pituitary gene expression

POMC was by far the most highly expressed gene in the pituitary, with an expression level an order of magnitude higher than the next most highly expressed gene. The robust expression of POMC was expected in anterior pituitary tissue, as it codes for a prohormone that is spliced into the primary hormonal product of both melanotrophs and corticotrophs (Drouin 2016). PCSK2, a convertase for processing POMC (Chrétien and Mbikay 2016), was also highly expressed. The third and fourth most highly expressed genes, EEF1A1, a translation elongation factor that serves to target actin mRNA to pseudopodial cell projections for translation (Liu ), and GNAS, an adenylyl cyclase regulator (Weinstein ), represented interesting findings in light of the importance of pseudopodium formation and cAMP regulation in differences between tame and aggressive fox pituitary function described by other analyses in this study. Of particular interest was the significant upregulation of EEF1A1 in the tame fox pituitary. Other highly expressed genes were related to transcription and mRNA splicing [DDX5/P68, SRRM2, DDX17, and HNRNPH1 (Jurica ; Lin ; Fuller-Pace and Ali 2008; Russo ; Wang )]. Along with the highly enriched GO terms associated with RNA processing and splicing for the top 0.1% most highly expressed genes, these findings suggest that the fox pituitary is highly active in RNA production. This enrichment is noteworthy due to the association of the Brown module with mRNA splicing functionality and the differentially spliced exons identified in the two populations of foxes.

Pseudopodia, cell motility, and networks

GO analysis identified pseudopodium organization as the most significant term associated with DEG in pituitary tissue, and an analysis separating upregulated genes in tame from aggressive samples identified enriched terms related to pseudopodium organization specifically in aggressive samples. In addition to the DEG included in this GO term, several genes known to be important for pseudopodium formation and functioning are not included in this term but are differentially expressed between pituitary samples of tame and aggressive foxes. These genes code for laminins (LAMA2 and LAMA5), a major component of the extracellular protein matrix, and integrins (ITGA1 and ITGA6), transmembrane receptors known to frequently interact with laminins to promote the formation of pseudopodia (Gu ), and to function in epithelial cell migration through pseudopod formation (Mercurio ). Two cadherin-family genes, CDH4 (R-cadherin or retinal cadherin) and CDH13 (T-cadherin, H-cadherin, or heart cadherin), are also differentially expressed; cadherins are implicated in cell motility, adhesion, and morphological changes, and these two cadherins particularly have been implicated in vascular remodeling (Dorrell ; Ivanov ). Weighted gene coexpression network analysis similarly described functionality related to cell migration and adhesion, particularly in the Antique White 4, Pale Turquoise, and Brown modules, with statistically significant increases in activity in the first two modules associated with aggressive samples. Significantly alternatively spliced genes in the tame vs. aggressive fox pituitary also demonstrated associations to integrins [MFAP5 (Gibson )] and cell adhesion [SORBS1 (Ribon ); EFNA5 (Davy and Robbins 2000); UNC5D (Jarjour )]. Together, these results demonstrate a recurrent theme of cell adhesion and pseudopodia formation throughout analyses of differential gene expression, GO term enrichment, weighted gene coexpression network analysis, and alternative splicing analysis, suggesting that this finding is a robust one.

Regulation of dendritic cell differentiation—another GO term highly significantly enriched for differentially expressed genes in the tame vs. aggressive pituitary—also contains genes associated with pseudopodia. Dendritic cells are antigen-presenting white blood cells that pass through the final stages of differentiation in tissues after being activated by exposure to antigen (Ziegler-Heitbrock ). Enrichment in genes related to regulation of their differentiation in the pituitary is a surprising finding; examination of the specific genes reveals two galectins related to cellular motility [LGALS3, LGALS9 (Sano ; Bi )] and an intracellular protein associated with cross-presentation of antigen [TMEM176B (Segovia )]. These findings may result from the migration of dendritic cells through the bloodstream to the pituitary gland and subsequent differentiation there. However, galectins and TMEM176B are not solely found in dendritic cells but are represented in numerous other tissues (Wada ; Stillman ; Cuajungco ). Their functionalities are not limited to dendritic cell differentiation; for example, galectin 3 has been associated with vesicle formation (Mehul and Hughes 1997) and localizes in dendritic cell lamellipodia (a type of pseudopodium) during migration (Hsu ); galectin 9 associates with integrins to regulate cell adhesion (Kasamatsu ; Nobumoto ). Therefore, while the enrichment of a GO term related to dendritic cell differentiation is difficult to explain in relationship to pituitary functionality, examination of specific DEG reveals functionality related to vesicle formation (associated with exocytosis) and pseudopodia formation and cell adhesion (associated with cellular networking of hormone-releasing cells).

Another highly enriched GO term, granulocyte chemotaxis, again describes a cell type (granulocytic white blood cell) that is not expected to be present in large quantities in the pituitary. This term contains two genes for receptors for interleukin 17A (IL-17), a pro-inflammatory cytokine, IL17RA and IL17RC (Gaffen 2009). While IL-17 is most commonly expressed in T cells, it is also widely expressed in other tissues, including epithelial tissue (Chesné ), which is the tissue type of the anterior pituitary (Conklin 1968). Interleukins signal through the JAK-STAT pathway (Subramaniam ; Sansone and Bromberg 2012), which is also represented by enriched GO terms (described below). IL-17 has different functions in different tissues, so its function in pituitary is difficult to ascertain, but it is known to enhance the production of tight junctions in epithelial tissue (Cua and Tato 2010). Tight junctions can control cell permeability and are associated with cell signaling and cell–cell adhesion (Steed ). In addition to receptors for the pro-inflammatory cytokine IL-17, this term also includes the pro- and anti-inflammatory cytokine TGFB2. Inflammatory cytokines, released during the acute phase of inflammation, have effects on the brain leading to behavioral changes. Specifically, they are known to cause “sickness behavior,” including social withdrawal (Parnet , Dantzer ). They also interact with the hypothalamic-pituitary-adrenal axis, both regulating, and being regulated by, glucocorticoid release (Parnet ). This term may therefore relate more to cell migration, adhesion, cytokine effects on behavior, and hormone release than to granulocyte chemotaxis in this tissue.

These findings are suggestive of differential function of pseudopodia, mediated by differences in integrin, laminin, and galectin expression, in tame and aggressive fox pituitaries. Pituitary lactotrophs and gonadotrophs are known to form pseudopodia. These cellular protrusions may function to both directly contact blood vessels and to increase cell motility to better position endocrine cells near blood vessels for hormone release (Navratil ). Pseudopodia may also function in the organization of anterior pituitary cellular networks for regulation of hormone release, allowing a small number of cells to coordinate the timing of hormone release to the bloodstream (Hodson ,b). These functions, while described in other pituitary endocrine cells, may be similar in ACTH-releasing corticotrophs. Notably, upregulation of pseudopodia-associated genes was largely observed in aggressive foxes, which are characterized by higher ACTH blood level relative to tame foxes.

Regulation of exocytosis

cAMP, a major mediator of exocytosis (Seino and Shibasaki 2005), was associated with a variety of DEG in the tame vs. aggressive fox pituitary. Of the 10 most significant GO terms associated with DEG in pituitary tissue, four are associated with regulation of cAMP. Other DEG not included in these GO groups are also associated with cAMP regulation, including two members of the family of phosphodiesterases involved in metabolism of cAMP: PDE4B and PDE7B. Interestingly, while PDE4B is upregulated in the tame pituitary, PDE7B is upregulated in the aggressive pituitary, suggesting a difference in function between the two proteins. Two members of the adenylyl cyclase family of proteins involved with synthesis of cAMP—ADCY7 and ADCY8—are both upregulated in the tame pituitary; along with phosphodiesterases, this family of proteins represents the major regulator of cAMP concentration in the cell. Adenylyl cyclase function is additionally regulated by the DEG CALCRL (a member of the Brown gene coexpression module), P2RY11 (also a member of the Brown module), and RAPGEF4 (Smith ; Qi ; Anand-Srivastava 2005; Sugawara ). A scaffolding protein, AKAP7 (a member of the Pale Turquoise module), is upregulated in the aggressive pituitary. AKAPs serve as anchoring proteins to spatially restrict cAMP activation (Kritzer ), and are also known to bind PDE4s to further spatially regulate cAMP (McCahill ; Dodge ; Taskén ). The cAMP regulation pathway therefore appears to be perturbed in tame pituitaries compared to aggressive pituitaries at several different levels. This supports the hypothesis that release of ACTH through cAMP-mediated exocytosis is an important factor in the different ACTH concentrations in tame and aggressive fox circulations.

The JAK-STAT pathway also appeared to be differentially regulated in tame vs. aggressive pituitaries; this pathway has multiple associations with cAMP and cell adhesion. Among the most significant GO terms (Table 1), four GO terms are related to the JAK-STAT pathway. The JAK-STAT pathway, besides its relationship to IL-17 as previously mentioned, is known to regulate cell migration in epithelial tissues (Pocha and Montell 2014). In the pituitary, JAK-STAT signaling operates synergistically with the cAMP pathway by mediating gp130 cytokine genes such as IL-6, LIF, IL-11, and CNTF to regulate POMC expression and ACTH release (Chesnokova and Melmed 2002) even in the absence of CRH (Kariagina ). Of the 11 DEG in the tame and aggressive fox pituitaries associated with this term, four (ERBB4, FLRT1, FLRT3, and LRRC4B) are also implicated in cell–cell adhesion and cellular migration (Gambarotta ; Jackson ; Yamagishi ; Wu ). All four are also upregulated in the aggressive pituitary. This pathway appears to link IL-17 signaling, associated with the highly significant enriched GO term granulocyte chemotaxis, with cAMP signaling.

Genes associated with secretory vesicles are also implicated in differences between tame and aggressive fox pituitaries. Two members of the synaptotagmin family of membrane trafficking proteins (Südhof 2002) are differentially expressed: SYT10, upregulated in the aggressive pituitary, and SYT13, upregulated in the tame pituitary. SYT10 is a calcium sensor required for exocytosis of secretory vesicles (Brose ); this gene is implicated in IGF-1 co-secretion with ACTH (Cao ), suggesting IGF-1 functionality in regulation of the stress response at the level of the pituitary (Eppler ). SYT10 is also a hub gene in the Antique White 4 module. SYT13 is an atypical synaptotagmin, lacking specific amino acid residues that are conserved in other family members, which may limit its calcium binding activity; though its function is not well understood, the SYT13 protein has been associated with constitutive vesicular transport, and shown to be upregulated in contextual fear conditioning (Han ). Another DEG associated with secretory vesicle function is VAMP1 (synaptobrevin), a membrane protein with binding specificity for vesicular proteins (Calakos ); this gene is upregulated in the tame pituitary. Overall, genes associated with exocytosis, both at the level of vesicle function and at the level of cAMP, display differences in expression between tame and aggressive fox pituitaries, not just individually, but in enriched GO terms and gene coexpression modules.

Regulation of ACTH release

The most significantly DEG, as well as the gene with the largest log2 fold change, was GRK7—a member of a gene family containing important regulators of CRH receptor function (Grammatopoulos 2012)—upregulated in the aggressive pituitary. Activation of the CRH receptor stimulates ACTH synthesis and release. Notably, however, the G protein-coupled receptor kinase (GRK) family is also associated with regulation of cell motility, including pseudopodium formation and cell adhesion (Chodniewicz and Klemke 2004; Penela ). GRK7 represents a potential candidate gene for future studies of the functional differences between the tame and aggressive fox pituitary.

The AVP gene, which codes for the vasopressin protein, is significantly upregulated in the tame fox pituitary. Vasopressin is known to potentiate ACTH release, but not POMC production, in the anterior pituitary (Levin ), both in endocrine fashion through release from the posterior pituitary into the systemic circulation, and in paracrine fashion through local synthesis and release (Denef 2008). However, interspecies analysis of adjacent, recently diverged genes such as AVP and OXT (Acher 1980) introduces the potential to conflate paralogy and homology, because sequence divergence has occurred between gene copies both within each species as well as within genes but between species. As a result, some fox OXT reads may have inappropriately mapped to dog AVP, resulting in apparent upregulation due to artifacts of the alignment algorithm rather than biological shifts. Therefore, the apparent upregulation of AVP should be treated cautiously until more work can more clearly define differences between dog and fox AVP sequences.

Protein targeting

Four highly significantly enriched GO terms (see Table 1) are related to protein targeting and signal recognition proteins (SRP), and genes in these terms appear to be upregulated in tame rather than aggressive pituitary tissues (see Table S10). These terms all describe processes controlling the translocation of proteins during translation to the organelle where they will ultimately reside. Notably, however, all but one of the DEG associated with these terms are associated with ribosomal proteins, i.e., with protein elongation rather than protein translocation. The remaining gene, ZFAND2B, is associated with protein folding as well as translocation (Glinka ). Interestingly, all of these genes are up regulated in the tame fox. This finding is surprising, as the aggressive pituitary is believed to release more ACTH and therefore might be expected to synthesize a greater number of proteins related to hormone synthesis and release.

Complexity of behavioral regulation

The observed differences between DEG in pituitaries from the tame and aggressive lines were less than twofold. A relatively small magnitude of changes was expected due to the outbred nature of these lines of foxes and complex genetic inheritance of their behavioral differences (Kukekova , Johnson ; Nelson ). Additionally, differences in gene expression between the two lines may not be a precise indicator of differences in protein expression, or of protein activity. A comparison of intraspecific differences in gene and protein between humans and chimpanzees suggested that protein expression levels are under stronger evolutionary constraints than are gene expression levels, and, therefore, a change in gene expression may not always indicate a corresponding change in protein expression (Khan ). Protein activity, moreover, may be modified not just by gene or protein expression changes, but by translational regulation or, post-translationally, of protein modifications such as phosphorylation, glycosylation, lipidation, ubiquitination, and acetylation. Such modifications may affect protein stabilization, activation, folding, and interaction with other proteins (Deribe ), and, thereby, may potentially significantly affect cellular processes without an accompanying perturbation of gene expression. Therefore, the transcriptome study presented here identified genes and gene pathways likely implicated in behavioral and physiological differences between tame and aggressive foxes, but further proteomic studies are needed to elucidate post-transcriptional differences between the two lines.

Conclusions

Experimentally domesticated foxes provide a novel model for the study of the behavioral phenotype of domesticated animals. An attenuated ACTH response to psychogenic stressors has been shown in the tame foxes, but the cellular mechanisms resulting in the inhibition of this response are not understood. Differences in ACTH concentrations in tame and aggressive fox blood have been demonstrated (Trut ), but both POMC mRNA and ACTH protein concentrations in tame and aggressive fox pituitaries are similar (Gulevich ). These results suggest that differences between the two lines may lie in regulation of ACTH release through exocytosis rather than regulation of ACTH synthesis and processing or in corticotroph cell number. In this study, examination of gene expression through sequencing of mRNA demonstrated differences in formation of pseudopodia and in signaling by cAMP. Pseudopodia formation have been implicated in release of pituitary hormones, such as lutenizing hormone, into blood vessels and in the formation of regulatory networks for coordinating hormone release by endocrine cells, while cAMP is a major regulator of exocytosis.

We have demonstrated here for the first time the possibility that genes in networks related to cAMP regulation, cell adhesion, and pseudopod formation may have been targeted in selection for the domesticated phenotype. Although our findings are in the experimentally domesticated fox, the similarity of phenotypes across domesticated species suggest the possibility of selection on a shared set of gene pathways. Therefore, our findings may provide new avenues for investigation of the biological differences between domesticated and wild species.

Supplementary Material

Supplemental material is available online at www.g3journal.org/lookup/suppl/doi:10.1534/g3.117.300508/-/DC1.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.