Genome-Wide Association Study on Three Behaviors Tested in an Open Field in Heterogeneous Stock Rats Identifies Multiple Loci Implicated in Psychiatric Disorders.

Many personality traits are influenced by genetic factors. Rodents models provide an efficient system for analyzing genetic contribution to these traits. Using 1,246 adolescent heterogeneous stock (HS) male and female rats, we conducted a genome-wide association study (GWAS) of behaviors measured in an open field, including locomotion, novel object interaction, and social interaction. We identified 30 genome-wide significant quantitative trait loci (QTL). Using multiple criteria, including the presence of high impact genomic variants and co-localization of cis-eQTL, we identified 17 candidate genes (Adarb2, Ankrd26, Cacna1c, Cacng4, Clock, Ctu2, Cyp26b1, Dnah9, Gda, Grxcr1, Eva1a, Fam114a1, Kcnj9, Mlf2, Rab27b, Sec11a, and Ube2h) for these traits. Many of these genes have been implicated by human GWAS of various psychiatric or drug abuse related traits. In addition, there are other candidate genes that likely represent novel findings that can be the catalyst for future molecular and genetic insights into human psychiatric diseases. Together, these findings provide strong support for the use of the HS population to study psychiatric disorders.

1. Introduction

Many personality traits are predictors of vulnerability to addiction (1). For example, individuals with symptoms of anxiety are more likely to be smokers (2, 3), and novelty seeking is positively correlated with both smoking onset (4) and cocaine abuse (5). In addition, the social environment plays a critical role in the development and treatment of addiction (6). Many of these phenomena can be modeled using rodents to unveil their neural, genetic, and molecular mechanisms (7–10).

The open-field test (OFT) is a widely used behavioral test for measuring anxiety-like and exploratory behavior in rodents (11–14). A rodent is typically placed in an open chamber surrounded by tall walls. Video recording of the rodent's locomotor movements is then analyzed. In general, rats spend most of the testing session walking along the wall (i.e., thigmotaxis). Increased time spent in the center of the area or decreased latency to enter the center are interpreted as indications of lower anxiety. The OFT is widely used to model anxiety and is sensitive to the anxiolytic-like effects of classical benzodiazepines, and 5-HT1A receptor agonists (11). The novel object interaction test (NOIT) is usually conducted in an open arena where a novel object is placed in the center. The time spent and distance traveled around the object zone are used as indicators of preference for novelty. Novel object interaction has been considered as an important predictor in addiction-like traits (15, 16) and high novelty preference increases the propensity for addictive drug-seeking behavior (9, 17, 18). There are multiple different methods for conducting social interaction test (SIT) in rats (19–21). In general, an unfamiliar stimulus rat and the rats to be tested are placed in the same arena. While manual scoring of social interaction often allows both rats to be freely moving, experiments using automated video analysis often limit the movement of the stimulus rat. Computer algorithm then extract the time spend and distance traveled by the test rat around the stimulus rat, which reflects the social tendency of the test rat.

The heterogeneous stock (HS) rats were originally derived from interbreeding eight inbred strains (22). An analysis on these founders reported 7.2 million single nucleotide variants (23). This population has been maintained as outbred for more than 90 generations. The chromosomes of individuals in this population represent a genetic mosaic of the founders' haplotypes, with the average distance between recombination events in the centiMorgan range (24). This allows for genetic mapping to only a few million bases (Mb), a much smaller region than what can be identified using traditional F2 intercross or backcross mapping strategies. Several high-resolution genome-wide association studies (GWAS) (23, 25–27) have been successfully carried out. Here we report the results on associations of genomic loci with measures obtained from OFT, NOIT and SIT. These analyses were based on an expanded data set that contained about twice the sample size of that reported previously (28). These data were collected as part of a larger GWAS on socially acquired nicotine intravenous self-administration, which will be the subject of a separate publication.

2. Materials and Methods

2.1. Animals

The N/NIH heterogeneous stock (HS) rat (RRID:RGD_2314009), was created at the NIH in 1984 by interbreeding the following eight inbred founder strains: ACI/N, BN/SsN, BUF/N, F344/N, M520/N, MR/N, WKY/N and WN/N (22). The HS rats used in this study were sent from The Medical College of Wisconsin to the University of Tennessee Health Science Center (UTHSC) at 3–6 weeks of age. A total of 16 batches of HS rats were transferred between October 27, 2014 and September 20, 2018. Each batch consisted of 25 males and 25 females that were used as breeders. After a 2-week quarantine period, rats were transferred to a reversed 12 h light-dark cycle (lights off at 9:00 a.m.) housing room. Breeding pairs were assigned according to an algorithm that maximized the genetic diversity of the offspring. Litters were culled to a maximum of 8 pups to ensure a consistent nutritional environment. Rats were weaned on postnatal day (PND) 21. A radio frequency identification (RFID) chip was inserted subcutaneously into each rat at the time of weaning. Two male and two female rats per litter were used for behavioral studies. Sprague-Dawley (SD) rats (20 for each sex, purchased from Harlan Laboratories, Madison, WI, RRID:RGD_737903) were used as the stimulus rats in the social interaction test. Teklad Irradiated LM-485 Mouse/Rat Diet and water were provided ad libitum. All rats were group-housed with 2–4 same-sex peers throughout the experiments to avoid social isolation. All procedures were conducted in accordance with the NIH Guidelines concerning the Care and Use of Laboratory Animals, as approved by the Institutional Animal Care and Use Committee of the University of Tennessee Health Science Center.

2.2. Study Design

All HS rats (626 males and 620 females in total from 16 batches) were adolescents when tests began. Their age was 31.8 ± 2.6 (mean ± STD) on the day of the OFT. Adolescent rats were used because the onset of many psychiatric diseases occur during this age (29). Each HS rat was tested in all three behavioral tests, one test per day, in the following sequence: OFT, NOIT, and SIT. All tests were conducted in the dark phase of the light cycle (9 a.m.–4 p.m.) and were conducted in the same open field and recorded using the same video capture system.

2.3. Behavioral Testing Procedure

2.3.1. Open Field Test

Two OFT arenas were constructed using black acrylic glass, measuring 100 cm (L) × 100 cm (W) × 50 cm (H), which were placed side by side. The floors were covered by wood boards painted with either black or white acrylic paint (ART-Alternatives, ASTM D-4236, Emeryville, CA, USA) to contrast the coat of the animals (i.e., a black board was used for rats with white fur). The test chambers were illuminated by a long-range, 850-nm infrared light (LIR850-70, LDP LLC, Carlstadt, NJ) located 160 cm above the center of the two test chambers. No source of visible light was present during behavioral testing, with the exception of a flat panel monitor (Dell 1908FP). A digital camera (Panasonic WV-BP334) fitted with an 830 nm infrared filter (X-Nite830-M37, LTP LLC, Carlstadt, NJ) and located next to the infrared light source was used to record the behavior of the rats. All rats were released at the same corner of the test chamber, and data were collected for 1 h.

2.3.2. Novel Object Interaction Test

This test was conducted the day after the OFT in the same arena. A cylindrical rat cage constructed using 24 aluminum rods (30 cm in length) spaced 1.7 cm apart was used as the novel object. The bottom and top of the cage (15 cm in diameter) were manufactured using a 3D printer from polylactic acid. The design can be downloaded from https://github.com/chen42/RatSocialInteractionTest. The novel object was placed in the center of the arena before testing. The test duration was 20 min and was recorded using the same camera as that used in the OFT.

2.3.3. Social Interaction Test

This test was conducted the day after the NOIT. This test compares the preference of a subject rat for a stimulus rat restricted in a cylindrical cage (i.e., the novel object used in the NOIT) against an empty cylindrical cage. The test arena was reduced to 100cm(L) × 60cm(W) × 50cm(H) by using a black board placed vertically in the arena. Two cylindrical cages described above were placed ~30 cm away from the walls on opposite sides (i.e., similar to the arrangement commonly used in the three-chamber test). A randomly selected stimulus Sprague-Dawley rat of the same sex and similar weight as the HS test rat was placed into one of the cylindrical cages (kept the same throughout the experiment) 5 min before the HS subject rat was placed into the arena. The stimulus and subject rats were never housed together and thus were unfamiliar to each other. No social isolation was conducted on either rat. Each stimulus rat was used no more than once per day. The test duration was 20 min and was recorded using the same camera as that used in the OFT.

2.3.4. Analysis of Video Data

Ethovision XT video tracking system (RRID:SCR_000441, Version 4.0, Noldus Information Technology, The Netherlands) was used to analyze the videos recorded in all behavioral tests. After identifying the arena and calibrating the size of the arena, specific zones in the arena were outlined. For OFT and NOIT, one center zone, which was a circular region with a diameter of 20 cm, was used. For the SIT, one object zone and one social zone, both were circular regions with diameters of 20 cm, corresponding to the two cylindrical cages, respectively, were specified. The extracted data included the total distance traveled in the arena, the duration and the frequency the test rat was present in specific zones, the distance of the subject to the zones, and the latency of the test rat entering the zones. The center of the subject rat was used for all calculations. Phenotypic correlations were determined using the Pearson test.

2.4. Pre-processing of Phenotype Data

All phenotype data were stored in the C-GORD (RRID:SCR_021866) relational database. For genetic analysis, each trait was quantile-normalized separately for males and females; this approach is similar to using sex as a covariate. Other relevant covariates (including age, batch number, and coat color) were identified for each trait, and covariate effects were regressed out if they were significant and if they explained more than 2% of the variance. Residuals were then quantile-normalized again, after which the data for each sex were pooled prior to further analysis. This approach removed mean differences due to sex; further, it did not attempt to model gene-by-sex interactions.

2.5. Genotyping and Estimates of Heritability

Genotypes were determined using genotyping-by-sequencing (GBS), as described previously (30). This produced approximately 3.5 million single nucleotide polymorphisms (SNP) with an estimated error rate <1%. Variants for X- and Y-chromosomes were not called. We used this set of SNPs for GWAS, genetic correlations, and heritability estimates. We used GCTA-GREML (31) analysis to estimate proportion of variance attributable to SNPs.

2.6. Genetic Mapping

GWAS analysis employed a linear mixed model, as implemented in the software GCTA (32), using a genetic relatedness matrix (GRM) to account for the complex family relationships within the HS population and the Leave One Chromosome Out (LOCO) method to avoid proximal contamination (33, 34). Significance thresholds were calculated using permutation. Because all traits were quantile normalized, we used the same threshold for all traits (35). To identify quantitative trait loci (QTL), we scanned each chromosome to determine if there was at least one SNP that exceeded the permutation-derived threshold of −log10(p) > 5.6, which was supported by a second SNP within 0.5 Mb that had a p-value that was within 2 − log10(p) units of the most significant SNP.

There could be more than one QTL on the same chromosome for one trait. We resolve the dependency and determine their locations as follows: we used the top SNP from the most significant QTL as a covariate and performed a second GWAS of the chromsome in question. If the resulting GWAS had an additional SNP with a p-value that exceeded our permutation-derived threshold, it was considered to be a second, independent locus. This process was repeated (including all previously significant SNPs as covariates), until no more QTLs were detected on a given chromosome. Linkage disequilibrium (LD) intervals for the identified QTL were determined by identifying those markers that had a high correlation coefficient with the peak marker (r2 = 0.6).

Genetic fine mapping were conducted using Credible Set analysis (36) and SuSieR (37). The analysis determines the 99% credible set by a Bayesian approach, that is the smallest set of SNPs in a genomic region that were 99% likely, to contain the causal SNPs. SuSieR also uses a Bayesian approach but also quantify uncertainty in which variants should be selected when multiple, highly correlated variants compete with one another.

We used fastENLoc (38) and a LD cutoff-based method to colocalize behavioral and gene expression QTLs. For the LD cutoff-based method, we retained those behavioral and gene expression QTLs where the top SNPs were in strong LD (i.e., r2 > 0.6). The gene expression data were collected from 88 naive adult HS rats. Five brain regions (prelimbic, infralimbic, and orbitofrontal cortex, lateral habenula, and nucleus accumbens core) were collected for RNA-seq from each rat (39).

3. Results

3.1. Sex Differences

We found that many of the traits measured in OFT, NOIT, and SIT are different between males and females (Supplementary Table S1). In OFT, with the exception of latency of entering the center zone, all traits have statistically significant sex differences. In addition, four out of six traits in NOIT and seven out of eleven traits in SIT are different between males and females. The range of effect size (Cohen's d) for statistically significant differences is (0.14, 0.31). Our genetic analysis quantile-normalized each trait separately for males and females. This approach removed mean differences due to sex and allowed us to combine males and females in the same analysis to increase the power of GWAS.

3.2. Phenotypic Correlations

We calculated Pearson correlation between the 23 traits (Figure 1). We found 197 correlations with un-adjusted p-values < 0.05. Most of these correlations have relatively low Person coefficient (mean is 0.23, median is 0.18). However, due to the large sample size, most of these correlations are highly significant (median −log10(p) is 7.8). In general, correlations of traits obtained from the same behavioral test are among the strongest. For example, frequency of visiting the center and duration of staying in the center are positively correlated in OFT (r = 0.76), and duration in the social zone and distance to the social zone in the SIT are negatively correlated (r = –0.76). Most of these correlations are expected from the definitions of these variables.

Figure 1

Heatmap showing the correlations between behavioral traits. The color scheme represents the direction of the correlation, whereas the intensity of the colors and the size of the circles are proportional to coefficients of the correlation. The cross signs indicates that the correlation of the two traits is not statistically significant (p > 0.05).

Among the correlations of variables derived from two different behavioral tests, correlations for measures of distance traveled are among the highest (range of Pearson r: 0.39–0.47, e.g., Figures 2A,B). Distance traveled in the OFT is also correlated with duration of center time in the NOIT (e.g., Figure 2C). Interestingly, the frequencies of visiting the center of the area in the NOIT is correlated with the frequency of visiting the social zone in the SIT (Figure 2D). In contrast, OFT center frequency is negatively correlated with NOIT mean distance to center in NOIT (Figure 2E), and distance to object zone in SIT is negatively correlated with center frequency in NOIT (Figure 2F).

Figure 2

Selected scatter plots for correlation between behavioral tests shown in Figure 1. (A) OFT vs. NOIT. (B) OFT vs. SIT. (C) OFT vs. NOIT. (D) SIT vs. NOIT. (E) OFT vs. NOIT. (F) SIT vs. NOIT.

3.3. Heritability

SNP heritability estimates (h2) for traits are provided in Table 1. In all the three behavioral tests, total travel distance has the highest heritability. In OFT, all heritability estimates are between 0.28 and 0.38, with the exception of that for latency of entering the center zone (h2 = 0.08). Heritability estimates for variables from the NOIT are slightly lower than that of the OFT; most of them are in the range of 0.21–0.29, with the exception of that for the latency of entering the center zone (h2 = 0.10). Heritability estimates for various measures of the SIT are in the range of 0.10–0.28. Interestingly, heritability estimates for measures on the social zone are consistently greater than those for the object zone.

Table 1

Heritability of open field (OFT), novel object (NOIT) and social interaction (SIT) tests.

Test	Trait	Heritability ±SE
OFT	Duration in center zone	0.284 ± 0.045
	Frequency of entering center zone	0.323 ± 0.044
	Latency of entering center zone	0.083 ± 0.034
	Mean distance to center zone	0.295 ± 0.043
	Total distance to center zone	0.300 ± 0.043
	Total travel distance	0.379 ± 0.044
NOIT	Duration in center zone	0.247 ± 0.043
	Frequency of entering center zone	0.209 ± 0.041
	Latency of entering center zone	0.100 ± 0.034
	Mean distance to center zone	0.249 ± 0.042
	Total distance to center zone	0.221 ± 0.041
	Total travel distance	0.287 ± 0.044
SIT	Duration in object zone	0.161 ± 0.037
	Duration in social zone	0.275 ± 0.040
	Frequency of entering object zone	0.177 ± 0.036
	Frequency of entering social zone	0.215 ± 0.036
	Latency of entering object zone	0.082 ± 0.032
	Latency of entering social zone	0.142 ± 0.034
	Mean distance to object zone	0.165 ± 0.038
	Mean distance to social zone	0.265 ± 0.041
	Total distance to object zone	0.153 ± 0.037
	Total distance to social zone	0.265 ± 0.041
	Total travel distance	0.281 ± 0.040

3.4. Identification of Multiple QTLs

In Table 2, we present SNPs that are significantly associated with the phenotypes. The genome-wide statistical significance of the association is determined by −log10P values greater than 5.609. For OFT, there are 9 significant loci for 5 traits. We did not find a significant QTL for Duration in center zone (h2 = 0.284±0.045). We identified two loci for Frequency of entering center zone and Total travel distance, 3 loci for Total distance to center zone. We found 4 NOIT traits have significant loci. Among them, Total distance to center zone has 3 loci and Mean distance to center zone has 2 loci. We did not find any significant loci for Frequency of entering center zone (h2 = 0.209 ± 0.041) and Latency of entering center zone (h2 = 0.100 ± 0.034). For SIT, we identified significant loci for all traits except Latency of entering object zone which has heritability of h2 = 0.082 ± 0.032. We found 2 loci for the traits Latency of entering social zone, Mean distance to social zone Total distance to social zone and Total travel distance. All genome-wide significant loci are shown in Figure 3. Genetic mapping of individual traits are shown as Manhattan plots as Supplementary Figures S1–S23. Regional association plots for representative traits are shown in Figures 4–6 for OFT, NOIT, and SIT, respectively. Other regional association plots are provided as Supplementary Figures S24–S50.

Table 2

QTL for open field (OFT), novel object interaction (NOIT), and social interaction (SIT) tests.

Test	Trait	Top SNP	−log10P	Interval size	Number of genes
OFT	Frequency of entering center zone	chr1:24043699	5.714	0.12 Mb	5
	Frequency of entering center zone	chr4:118013062	5.777	2.0 Mb	47
	Latency of entering center zone	chr8:120910798	5.609	1.0 Mb	12
	Mean distance to center zone	chr4:58009499	7.469	2.4 Mb	60
	Total distance to center zone	chr4:58009499	7.254	2.4 Mb	60
	Total distance to center zone	chr4:118013062	6.099	2.0 Mb	47
	Total distance to center zone	chr14:44904830	5.741	2.1 Mb	44
	Total travel distance	chr10:94549701	7.286	4.2 Mb	98
	Total travel distance	chr11:33359859	8.268	0.92 Mb	23
NOIT	Duration in center zone	chr4:112234344	6.028	1.2 Mb	8
	Mean distance to center zone	chr4:112234344	6.598	1.2 Mb	8
	Mean distance to center zone	chr6:119975012	5.692	0.95 Mb	3
	Total distance to center zone	chr1:144080083	5.969	4.1 Mb	109
	Total distance to center zone	chr4:112234344	5.975	1.2 Mb	8
	Total distance to center zone	chr4:156801420	5.622	4.4 Mb	127
	Total travel distance	chr6:120117521	5.640	0.95 Mb	3
SIT	Duration in object zone	chr18:65869186	6.414	3.4 Mb	22
	Duration in social zone	chr4:151128675	5.820	2.9 Mb	34
	Frequency of entering object zone	chr13:90335374	5.827	1.1 Mb	46
	Frequency of entering social zone	chr1:239076581	7.273	0.27 Mb	6
	Latency of entering social zone	chr10:52831274	6.052	0.34 Mb	2
	Latency of entering social zone	chr17:58611795	6.104	0.86 Mb	1
	Mean distance to object zone	chr19:20666789	6.746	1.6 Mb	23
	Mean distance to social zone	chr19:55339863	6.661	0.68 Mb	16
	Mean distance to social zone	chr4:150582701	5.884	1.1 Mb	19
	Total distance to object zone	chr19:20667417	6.619	1.6 Mb	23
	Total distance to social zone	chr19:55339863	6.643	0.68 Mb	16
	Total distance to social zone	chr4:150582701	5.788	1.1 Mb	19
	Total travel distance	chr14:34908176	5.648	0.74 Mb	10
	Total travel distance	chr14:41727329	5.627	0.85 Mb	5

Figure 3

Association of approximately 3 million SNPs with behavioral traits measured in OFT, NOIT, or SIT. The red horizontal line denotes the p-value for reaching genome-wide significance. The downward arrows denote the SNPs with the largest –log10(P) for each genome-wide significant association.

Figure 4

Regional association plot for frequency of entering center zone in OFT at chr1:24043699.

Figure 6

Regional association plot for latency of entering social zone in SIT at chr17:58611795.

Regional association plot for mean distance to center zone in NOIT at chr6:119975012.

3.5. Pleiotropic Loci

To determine if traits that mapped to the same location are pleiotropic, we considered the minor allele frequency (MAF), and the strain distribution pattern (SDP) of the most significant SNP among the 8 founder strains that were used to create the HS. Using these criteria, we did not observe any pleiotropic loci between the traits analyzed in different tests. However, we did identify pleiotropic loci between the traits of the same behavior test. Most of these traits are highly correlated, as shown in Figure 1. With the exception of three sets of QTL (Supplementary Table S3), all others share the same top SNP (Supplementary Table S2).

3.6. Candidate Gene Identification

The number of genes within the identified QTL ranges from 1 to 127 (mean: 30.1, median: 19, Table 2). There is only one region that contains a single gene: Adarb2 within chr17:58Mb for latency of entering social zone in SIT. However, it is also possible that the causal allele is a regulatory variant that is located in this interval but regulates a gene outside of the identified interval.

All other loci contained more than one gene. To identify candidate genes, we combined several criteria: (1) located in the credible set identified by either one of the fine mapping methods. (2) the presence of moderate or high impact variants located within the gene, as predicted by SnpEff (40). We also require these variants are in high LD with the top SNP. We identified 149 coding variants within 30 QTL, 8 of which were predicted to have a high impact (Supplementary Table S4). (3) the presence of a significant cis-eQTL in one or more of the five brain regions in a dataset containing 88 navie adult HS rats (39), (4) has a human ortholog that has been reported to be associated with psychiatric diseases (including drug abuse). When multiple candidates are present using the above criteria, we remove the gene with very low expression levels across all five regions in the RNA-seq data set (e.g., FPKM < 0.5) and select the candidate with the strongest support for the literature. Combining these criteria with a literature search conducted using GeneCup (41), we identified plausible candidate genes within 17 loci (Table 3).

Table 3

Candidate genes.

Test	Trait	Top SNP	Candidate gene	Supporting evidence	Human GWAS	Expression level (FPKM)	Gene function
OFT	Frequency of entering center zone, Total distance to center zone	chr4:118013062	Cyp26b1	Missense variants, cis-eQTL in IL and PL	Schizophrenia (42, 43)	IL 7.29 ± 2.26	Inactivate all-trans retinoic acid (44)
OFT	Total distance to center zone	chr14:44904830	Fam114a1	Missense variants, cis-eQTL in LHb	Alcohol consumption measurement (45)	LHb 3.72 ± 0.77	Also known as Noxp20, neuronal cell development (46)
OFT	Total travel distance	chr10:94549701	Cacng4	cis-eQTL in Acbc	Bipolar disorder and Schizophrenia (47)	Acbc 55.22 ± 4.79	Calcium channel (48)
OFT	Total distance to center zone	chr4:58009499	Ube2h	cis-eQTL in OFC	Unipolar depression, mood disorder (49)	OFC 49.22 ± 2.72	Ubiquination of proteins (50)
NOIT	Duration in center zone, distance to center zone	chr4:112234344	Eva1a	Missense variants, cis-eQTL in PL, IL and OFC		LHb 6.74 ± 1.69	Formation of the autophagosome (51)
NOIT	Total distance to center zone	chr1:144080083	Sec11a	cis-eQTL in PL and IL	Unipolar depression, depressive symptom measurement, response to ketamine, bipolar disorder, schizophrenia (43, 52, 53)	Acbc 33.34 ± 3.73	Metabolism of proteins (54)
NOIT	Total distance to center zone	chr4:156801420	Mlf2	cis-eQTL in PL	Smoking status measurement (55)	Acbc 307.85 ± 40.01	Molecular chaperone in multi-protein complex assembly, signaling transduction, and endocytosis (56).
SIT	Duration in object zone	chr18:65869186	Rab27b	missense variants, cis-eQTL in LHb	Unipolar depression, bipolar disorder (57–59)	Acbc 8.03 ± 2.7	Vesicular fusion and trafficking (60)
SIT	Duration in social zone, distance to social zone	chr4:151128675	Ankrd26	Missense variants	smoking initiation (61)	IL 4.75 ± 1.21	Cell signaling (62)
SIT	Frequency of entering object zone	chr13:90335374	Kcnj9	Missense variants, cis-eQTL in IL, PL and OFC	Alcohol consumption measurement (63)	OFC 65.46 ± 6.64	Adult neurogenesis (64), cocaine addiction (65)
SIT	Frequency of entering object zone	chr1:239076581	Gda	cis-eQTL in IL	General cognitive function (66)	Acbc 84.67 ± 16.92	Cypin, cytoplasmic PSD95 Interactor (67)
SIT	Latency of entering social zone	chr10:52831274	Dnah9	Missense variants	Schizophrenia (68)	LH 3.59 ± 1.33	Component of microtubule (69)
SIT	Latency of entering social zone	chr17:58611795	Adarb2	cis-eQTL in Acbc and LH	Unipolar depression, smoking status measurement, systolic blood pressure (58, 70)	PL 2.22 ± 0.68	Editing of neurotrasmiter mRNA (71)
SIT	Distance to social zone	chr19:55339863	Ctu2	Missense variants, cis-eQTL in IL	Autism spectrum disorder symptom (72)	LHb 6.39 ± 1.14	Post-transcriptional modification of tRNAs (73)
SIT	Distance to social zone	chr4:150582701	Cacna1c		Schizophrenia (74), biopolor disorder (75)	PL 6.76 ± 1.49	Calcium channel (76)
SIT	Total travel distance	chr14:34908176	Clock	Missense variants, cis-eQTL in LHb		Acbc 11.15 ± 1.82	Regulate circadian rhythms (77, 78)
SIT	Total travel distance	chr14:41727329	Grxcr1	cis-eQTL in Acbc and LHb	Cognitive decline in depression (79)	IL 2.97 ± 1.02	S-glutathionylation of proteins (80)

Acbc, Accubens core; IL, infralimbic cortex; LHb, lateral habenular; OFC, orbitofrontal; PL, prelimbic cortex.

In addition, for total distance to the novel object zone, the QTL region on chr1 (144 Mb, size: 4.1 Mb, Supplementary Figure S4) contains 69 gene with human orthologs. We found 14 of these genes have been reported in human GWAS to be associated with psychiatric conditions or addiction with genome-wide significance (ACAN, ADAMTSL3, ALPK3, CPEB1, FES, FURIN, LINC00933, MIR9-3HG, MRPL46, NMB, POLG-DT, SEC11A, ZNF592, and ZSCAN2, Supplementary Table S5). Three additional genes with sub-threshold significance in human GWAS are also included. These genes are all located in a syntenic region on human chromosome 15 (82.5–90.8 Mb). Although based on the criteria described above, Sec11a is the best candidate gene (Table 3), it is possible that this region contains multiple genes that are associated with the trait.

4. Discussion

As part of a GWAS on intravenous nicotine self-administration in adolescent HS rats that we are conducting (28, 81), we collected several behavioral phenotypes related to anxiety, novelty exploration, and social interaction. We have previously reported that these behavioral traits contribute to the variation in nicotine intake (28). We report here GWAS results of three behavioral traits: OFT, NOIT, and SIT, which were all conducted in the same open field. We identified 30 QTLs for 23 traits. Using a set criteria outlined above, we identified 17 candidate genes.

OFT, NOIT, and SIT are widely used behavioral assays in rodents. With over 1,200 rats, ours represent some of the largest data collected using these assays. Similar to our interim report on this data set (28), we found a large number of correlations with relatively low coefficients (e.g., r < 0.4) but with high statistical significance. It is likely that these correlated traits are controlled by the same behavioral processes and thus are influenced by the same genetic factors. In fact, our genetic analysis did find several pleiotropic sites (Supplementary Table S3). Almost all pleiotropic loci are reported for traits measured in the same behavior assay. It is likely that further increasing sample size will provide greater statistical power to detect pleiotropic effect across different behavioral assays.

Many of the candidate genes in this study have been associated with psychiatric or drug abuse traits in humans. For example, we identified Cyp26b1, a retinoic acid degrading enzyme, as a candidate gene for the frequency of entering the center zone and total distance to the center zone in OFT; both of which are measures of anxiety-like behaviors (rats with more anxiety-like behavior would enter the center zone less frequntly and have smaller distance to the center zone) (11). Cyp26b1 has been associated with Schizophrenia in several human GWAS (42, 43). Anxiety symptoms are common in schizophrenia patients (82, 83). Cyp26b1 is expressed in parvalbumin-positive interneurons (84). Most interestingly, knockdown Cyp26b1 in the nucleus accumbens shell decreased anxiety-like behavior (44).

Among the candidate genes for NOIT, Eva1a is a candidate gene for the duration stayed in the center zone that contained the novel object. Eva1a is supported by strong cis-eQTL and a missense variant but has no literature support. Thus, further evaluating the role of Evala could potentially lead to new mechanisms for novelty seeking-like behavior. Sec11a, a candidate gene for total distance in the center zone, is associated with depression and schizophrenia (43, 52, 53). Mlf2, a candidate gene for total distance to center zone in NOIT, is associated with smoking in humans (55) and has very high expression levels in the accumbens (Table 3).

For the SIT, we identified Cacna1c, encoding the Ca1.2 subunit of the L-type Ca2+ channel, as a candidate gene for distance to the social zone, where the stimulus rat resided. Cacna1c has been associated with schizophrenia (74) and bipolar disorder (75) in human GWAS. Both schizophrenia and bipolar disorders are associated with impairments in a range of social deficits (85, 86). In animal studies, Sprague-Dawley rats with heterozygotic deletion of the Cacna1c gene (homozygotic mutation is lethal) showed many deficits in social behavior. These included reduced levels of ultrasonic vocalizations during rough-and-tumble play, as well as social approach behavior elicited by playback of ultrasonic vocalizations (87, 88). In mice, a knockdown of Cacna1c in the nucleus accumbens significantly increased susceptibility to social stress (89). Knocking down of Cacna1c in the prefrontal cortex of adult mice also recapitulated many of the social deficits (90). Importantly, some of the behavioral effects of Cacna1c appear to interact with genetic background (91).

Among the other candidate genes for the SIT traits, Rab27b is involved in the presynaptic mechanism of long-term potentiation (92) as well as myelin biogenesis in oligodendrocytes (93). Ankrd26 is expressed in the arcuate and ventromedial nuclei and in the ependymaGda, also known as Cypin, in located in the postsynaptic density (95). Ctu2 is involved in post-translational modification of tRNAs (73). Adarb2 has been associated with home cage activity (96) and unipolar depression (58). The Clock gene is involved in the maintenance of locomotor rhythms (97). Mutations of the CLOCK gene have been implicated in many psychiatric disorders (98). Although these candidates are well supported by multiple lines of evidence, additional work is needed to confirm their causal relationship to the corresponding behavioral traits.

The total distance to the novel object zone is associated with chr1:144080083 (allele frequency: 0.91, −log10(p) = 5.969, size of interval: 4.1 Mb, Supplementary Figure S34). This SNP is also associated with the duration rats stayed in the novel object zone, although the p-value did not reach genome-wide significance (-logP = 4.63). This region contains 69 known genes. Its syntenic region on human Chr15 (82.5-90.8 Mb) is a hotspot for human pyschiatric diseases, containing 30 SNPs and 14 genes (ACAN, ADAMTSL3, ALPK3, CPEB1, FES, FURIN, LINC00933, MIR9-3HG, MRPL46, NMB, POLG-DT, SEC11A, ZNF592, and ZSCAN2) associated with generalized anxiety disorder, schizophrenia, bipolar disorder, obsessive compulsive disorder, attentions deficit hyperactivity disorder, autism spectrum disorder, and unipolar depression, smoking behavior, etc. These results are reported in 21 publications (Supplementary Table S5). Using the criteria described above, we identified Sec11a as the best candidate gene (Table 3). However, given the large number of genetic variants reported in human GWAS that are associated with psychiatric conditions within this syntenic region, it is very likely that this region contains multiple genes that are associated with novelty seeking-like behavior.

We include overlapping with human psychiatric GWAS results as part of the criteria in prioritizing candidate genes. It is possible that this approach could introduce bias and prevent us from making novel discoveries. For example, two (Cyp26b1 and Cacng4) of the four candidate genes for OFT have been associated with schizophrenia, rather than anxiety. However, many genetic variants are pleiotropic for multiple psychiatric diseases (99). For example, polygenic risk scores for schizophrenia have been associated with many other psychiatric diseases, such as anxiety disorder (100) or major depressive disorder (101), or cognitive performance (102). Together with other evidence, we believe considering human psychiatric GWAS findings when identifying candidate genes in our study, even when the behavior trait in rats does not map directly to the psychiatric disease, is still valid and will likely increase the translational value of our findings.

The presence of cis-eQTL in the brain is one of the strongest pieces of evidence that we use to prioritizes candidate genes. Fourteen of the 17 candidate genes we identified have cis-eQTL. Seven of these genes also contained missense mutations, which further support their biological function related to the phenotype they are associated with. The only candidate gene that is not support by either cis-eQTL or missense mutation is Cacna1c, associated with distance to the social zone. However, the role of Cacna1c in social behavior has been well documented in rats and mice (87–90). We also required all candidate genes located in regions confirmed by fine mapping. For example, the Crhr1 gene, which encodes corticotrophin release hormone receptor 1, is located in the locus for total travel distance in the OFT but is not supported by fine mapping. Dispite strong literature support for the role of Crhr1 in anxiety-like behavior in OFT in rats (103, 104) and mice (105), we nominated a different gene, Cacng4, to be the candidate gene for this locus. We anticipate further improvement in the statistical power of eQTL data and the availability of additional functional genomics data, such as 3D chromatin interaction (106, 107), will help us to identify additional candidate genes.

The HS rat population has already been successfully used in genetic mapping studies of physiological or behavioral traits (24, 25, 27). Prior study mapped several anxiety-like traits using zero maze (23). GWAS using HS to study behavioral regulation (108), response to cocaine cues (109), cocaine (110), nicotine (28, 81), or oxycodone (111) self-administration are underway. Our study adds to the literature 30 QTLs and 17 candidate genes for psychiatric related behavioral traits. Although we prioritized candidate gene selection based on functional genomics evidence, most of the candidate genes we identified have strong literature support for their role in human psychiatric diseases. This suggests that the rest of the candidate genes likely represent novel findings that can be the catalyst for future molecular and genetic insights on psychiatric diseases. In addition, these findings provide strong support for the use of the HS population in study psychiatric disorders.

Data Availability Statement

The original contributions presented in the study are publicly available. This data can be found through the C-GORD database at doi: 10.48810/P44W2 and through https://www.genenetwork.org.

Ethics Statement

The animal study was reviewed and approved by University Tennessee Health Science Center IACUC.

Author Contributions

HC and AP designed the study. TW and AG collected the data. AC, OP, and MG analyzed the data. MG, AP, and HC wrote the manuscript. All authors contributed to the article and approved the submitted version.

Funding

This work was supported by the National Institute on Drug Abuse (P50 DA037844).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.