A Genomic Region Containing REC8 and RNF212B Is Associated with Individual Recombination Rate Variation in a Wild Population of Red Deer (Cervus elaphus)

Recombination is a fundamental feature of sexual reproduction, ensuring proper disjunction, preventing mutation accumulation and generating new allelic combinations upon which selection can act. However it is also mutagenic, and breaks up favorable allelic combinations previously built up by selection. Identifying the genetic drivers of recombination rate variation is a key step in understanding the causes and consequences of this variation, how loci associated with recombination are evolving and how they affect the potential of a population to respond to selection. However, to date, few studies have examined the genetic architecture of recombination rate variation in natural populations. Here, we use pedigree data from ∼ 2,600 individuals genotyped at ∼ 38,000 SNPs to investigate the genetic architecture of individual autosomal recombination rate in a wild population of red deer (Cervus elaphus). Female red deer exhibited a higher mean and phenotypic variance in autosomal crossover counts (ACC). Animal models fitting genomic relatedness matrices showed that ACC was heritable in females (h2 = 0.12) but not in males. A regional heritability mapping approach showed that almost all heritable variation in female ACC was explained by a genomic region on deer linkage group 12 containing the candidate loci REC8 and RNF212B, with an additional region on linkage group 32 containing TOP2B approaching genome-wide significance. The REC8/RNF212B region and its paralogue RNF212 have been associated with recombination in cattle, mice, humans and sheep. Our findings suggest that mammalian recombination rates have a relatively conserved genetic architecture in both domesticated and wild systems, and provide a foundation for understanding the association between recombination loci and individual fitness within this population.

Genomic studies in humans, cattle, sheep and mice have shown that variation in recombination rate is often heritable, and may have a conserved genetic architecture (Kong et al. 2014;Ma et al. 2015;Johnston et al. 2016;Petit et al. 2017). The loci RNF212, REC8 and HEI10, among others, have been identified as candidates driving variation in rate, with PRDM9 driving recombination hotspot positioning in mammals (Baudat et al. 2010;Baker et al. 2017). This oligogenic architecture suggests that recombination rates and landscapes have the potential to evolve rapidly under different selective scenarios, in turn affecting the rate at which populations respond to selection (Barton and Charlesworth, 1998;Burt, 2000;Otto and Barton, 2001;Gonen et al. 2017). However, it remains unclear how representative the above studies are of recombination rate variation and its genetic architecture in natural populations. For example, experimental and domesticated populations tend to be subject to strong selection and have small effective population sizes, both of which have been shown theoretically to indirectly select for increased recombination rates to escape Hill-Robertson interference (Otto and Barton 2001;Otto and Lenormand 2002; but see Muñoz-Fuentes et al. 2015). Therefore, it may be that prolonged artificial selection results in different recombination dynamics and underlying genetic architectures. As broad recombination patterns are characterized in greater numbers of natural systems Theodosiou et al. 2016;Kawakami et al. 2017), it is clear that broad and fine-scale recombination rates and landscapes can vary to a large degree even within closely related taxa (Stapley et al. 2017). Therefore, determining the genetic architecture of recombination rate in nonmodel, natural systems are key to elucidating the broad evolutionary drivers of recombination rate variation and quantifying its costs and benefits at the level of the individual.
In this study, we investigate the genetic basis of recombination rate variation in a wild population of red deer (Cervus elaphus) on the island of Rùm, Scotland (Clutton-Brock et al. 1982). This population has been subject to a long term study since the early 1970s, with extensive pedigree and genotype information for 2,600 individuals at .38,000 SNPs (Huisman et al. 2016;Johnston et al. 2017). We use this dataset to identify autosomal crossover rates and their genetic architecture in .1,300 individuals. The aims of the study are to: (a) determine which common environmental and individual effects, such as age, sex and birth year affect individual recombination rates; (b) determine if recombination rate is heritable; and (c) identify genomic regions that are associated with recombination rate variation. Addressing these objectives will provide a foundation for future studies investigating the association between the genetic architecture of recombination rate and individual fitness, to determine how this trait evolves within contemporary natural populations.

Study population and genomic dataset
The study population of red deer is situated in the North Block of the Isle of Rùm, Scotland (57°02'N, 6°20'W) and has been subject to individual monitoring since 1971 (Clutton-Brock et al. 1982). Research was conducted following approval of the University of Edinburgh's Animal Welfare and Ethical Review Body and under appropriate UK Home Office licenses. DNA was extracted from neonatal ear punches, cast antlers and post-mortem tissue (see Huisman et al. 2016 for full details). DNA samples from 2880 individuals were genotyped at 50,541 SNP loci on the Cervine Illumina BeadChip (Brauning et al. 2015) using an Illumina genotyping platform (Illumina Inc., San Diego, CA, USA). SNP genotypes were scored using Illumina GenomeStudio software, and quality control was carried out using the check.marker function in GenABEL v1.8-0 (Aulchenko et al. 2007) in R v3.3.2, with the following thresholds: SNP genotyping success .0.99, SNP minor allele frequency .0.01, and ID genotyping success .0.99, with 38,541 SNPs and 2,631 IDs were retained. There were 126 pseudoautosomal SNPs identified on the X chromosome (i.e., markers showing autosomal inheritance patterns; Johnston et al. 2017). Heterozygous genotypes within males at non-pseudoautosomal X-linked SNPs were scored as missing. A pedigree of 4,515 individuals has been constructed using microsatellite and SNP data using the software Sequoia (see Huisman, 2017). The genomic inbreeding coefficient (F III ), was calculated for each deer in the software GCTA v1.24.3 (Yang et al. 2011), using information for all autosomal SNP loci passing quality control. A linkage map of 38,083 SNPs has previously been constructed, with marker orders and estimated base-pair positions known for all 33 autosomes (CEL1 to CEL33) and the X chromosome (CEL34)  and data archive doi: 10.6084/m9.figshare.5002562). All chromosomes are acrocentric with the exception of one metacentric autosome (CEL5).

Quantification of meiotic crossovers
A standardized sub-pedigree approach was used to identify the positions of meiotic crossovers . The full pedigree was split as follows: for each focal individual (FID) and offspring pair, a subpedigree was constructed that included the FID, its mate, parents and offspring ( Figure S1), where all five individuals were genotyped on the SNP chip. This pedigree structure allows phasing of SNPs within the FID, characterizing the crossovers occurring in the gamete transferred from the FID to the offspring. All remaining analyses outlined in this section were conducted in the software CRI-MAP v2.504a (Green et al. 1990) within the R package crimaptools v0.1  implemented in R v3.3.2. Mendelian incompatibilities within sub-pedigrees were identified using the prepare function and removed from all affected individuals; sub-pedigrees containing more than 0.1% mismatching loci between parents and offspring were discarded. The chrompic function was used to identify the grand-parental phase of SNP alleles on chromosomes transmitted from the FID to the offspring, and to provide a sex-averaged linkage map. Switches in phase indicated the position of a crossover ( Figure S1). Individuals with high numbers of crossovers per gamete (.60) were assumed to have widespread phasing errors and were removed from the analysis; the maximum number of crossovers for an individual was 45 in the remaining dataset.
Errors in determining allelic phase can lead to incorrect calling of double crossovers (i.e., $ 2 crossovers occurring on the same chromosome) over short map distances. To reduce the likelihood of calling false double crossover events, phased runs consisting of a single SNP were recoded as missing (390 out of 7652 double crossovers; Figure S2) and chrompic was rerun. Of the remaining double crossovers, those occurring over distances of # 10cM (as measured by the distance between markers immediately flanking the double crossover) were recoded as missing (170 out of 6959 double crossovers; Figure S2). After this process, 1341 sub-pedigrees passed quality control, characterizing crossovers in gametes transmitted to 482 offspring from 81 unique males and 859 offspring from 256 unique females.

Genetic architecture of recombination rate variation
Heritability and cross-sex genetic correlation: Autosomal crossover count (ACC) was modeled as a trait of the FID. A restricted maximum-likelihood (REML) "animal model" approach (Henderson, 1975) was used to partition phenotypic variance and examine the effect of fixed effects on ACC; these were implemented in ASReml-R (Butler et al. 2009) in R v3.3.2. The additive genetic variance was calculated by fitting a genomic relatedness matrix (GRM) constructed for all autosomal markers in GCTA v1.24.3 (Yang et al. 2011); the GRM was adjusted assuming similar frequency spectra of genotyped and causal loci using the argument -grm-adj 0. There was no pruning of related individuals from the GRM (i.e., we did not use the -grm-cutoff argument) as there is substantial relatedness within the population, and initial models included parental effects and common environment which controls for effects of shared environments between relatives. ACC was modeled first using a univariate model: where y is a vector of ACC; X is an incidence matrix relating individual measures to a vector of fixed effects, b; Z 1 ; and Z r are incidence matrices relating individual measures with additive genetic and random effects, respectively; a and u r are vectors of GRM additive genetic and additional random effects, respectively; and e is a vector of residual effects. The narrow-sense heritability h 2 was calculated as the ratio of the additive genetic variance to the sum of variance components estimated for all random effects. Model structures were tested with several fixed effects, including sex,F III and FID age; random effects included individual identity (i.e., permanent environment) to account for repeated measures in the same FID, maternal and paternal identity, and common environment effects of FID birth year and offspring birth year. A dominance GRM constructed in GCTA was also fitted but was not significant and explained ,1% of the phenotypic variance, and was therefore not included in the final models. The significance of fixed effects was tested with a Wald test, and the significance of random effects was calculated using likelihood-ratio tests (LRT, defined as 2 · the difference in log-likelihoods between the two models, distributed as x 2 with 1 degree of freedom) between models with and without the focal random effect.F III and individual identity were retained in all models, irrespective of statistical significance, to account for possible underestimation of ACC and pseudoreplication, respectively. As the variance in recombination rates differed between the sexes, models were also run for each sex separately.
Bivariate models of male and female ACC were run to determine whether additive genetic variation was associated with sex-specific variation and the degree to which this was correlated between the sexes. The additive genetic correlation r A was determined using the CORGH error-structure function in ASReml-R (correlation with heterogeneous variances) with r A set to be unconstrained. Model structure was otherwise the same as for univariate models. To determine whether genetic correlations were significantly different from 0 and 1, the unconstrained model was compared with models where r A was fixed at values of 0 or 0.999. Differences in additive genetic variance in males and females were tested by constraining both to be equal values using the CORGV error-structure function in ASReml-R. Models then were compared using LRTs with 1 degree of freedom.
Genome-wide association study: Genome-wide association studies (GWAS) of ACC were conducted using the function rGLS in the R library RepeatABEL v1.1 (Rönnegård et al. 2016) implemented in R v3.3.2. This function accounts for population structure by fitting the GRM as a random effect, and allows fitting of repeated phenotypic measures per individual. Models were run including sex andF III as fixed effects; sex-specific models were also run. Association statistics were corrected for inflation due to population stratification that was not captured by the GRM, by dividing them by the genomic control parameter l, which was calculated as the observed median x 2 statistic divided by the null expectation median x 2 statistic (Devlin et al. 1999). The significance threshold after multiple testing was calculated using a linkage disequilibrium (LD) based approach in the software K effective (Moskvina and Schmidt, 2008) specifying a sliding window of 50 SNPs. The effective number of tests was calculated as 35,264, corresponding to a P value of 1:42 · 10 206 at a = 0.05. GWAS of ACC included the X chromosome and 458 SNP markers of unknown position. It is possible that some SNPs may show an association with ACC if they are in LD with polymorphic recombination hotspots (i.e., associations in cis), rather than SNPs associated with recombination rate globally across the genome (i.e., associations in trans). Therefore, we repeated the GWAS modeling trans variation only, by examining associations between each SNP and ACC, minus the crossovers that occurred on the same chromosome as the SNP. For example, if the focal SNP occurred on linkage group 1, association was tested with ACC summed over linkage groups 2-33. Marker positions are known relative to the cattle genome vBTA_vUMD_3.1; in cases of significant associations with recombination rate, gene annotations and positions were obtained from Ensembl (Cattle gene build ID BTA_vUMD_3.1.89). LD was calculated between loci in significantly associated regions using the allelic correlation r 2 in the R package LDheatmap v0.99-2 (Shin et al. 2006) in R v3.3.2.
Regional heritability analysis: As a single locus approach, GWAS has reduced power to detect variants with small effect sizes and/or low linkage disequilibrium with causal mutations (Yang et al. 2011). Partitioning additive genetic variance within specific genomic regions (i.e., a regional heritability approach) incorporates haplotype effects and determines the proportion of phenotypic variance explained by defined regions. The additive genetic variance was partitioned across all autosomes in sliding windows of 20 SNPs (with an overlap of 10 SNPs) as follows (Nagamine et al. 2012;Bérénos et al. 2015): where y is a vector of ACC; X is an incidence matrix relating individual measures to a vector of fixed effects, b; v i is a vector of additive genetic effects explained by autosomal genomic region in window i; v ni is the vector of the additive genetic effects explained by all autosomal markers not in window i; Z 1 ; and Z 2 are incidence matrices relating individual measures with additive genetic effects for the focal window and the rest of the genome, respectively; Z r is an incidence matrix relating individual measures with additional random effects, where u r is a vector of additional random effects; and e is a vector of residual effects. The mean window size was 1.29 6 0.32 Mb. Models were implemented in ASReml-R (Butler et al. 2009) in R v3.3.2. GRMs were constructed in the software GCTA v1.24.3 with the argument -grm-adj 0 (Yang et al. 2011). The significance of additive genetic variance for window i was tested by comparing models with and without the Z 1 v i term with LRT (x 2 1 ). We attempted to model the X-chromosome GRMs as calculated in GCTA omitting the -grm-adj 0 argument. However, X-chromosome models had negative pivots in the GRMs, possibly a result of small sample and window sizes, and/or due to the fact that we could not modify the assumed frequency spectra of causal and typed loci as we could for the autosomal markers. To correct for multiple testing, a Bonferroni approach was used, taking the number of windows and dividing by 2 to account for window overlap; the threshold P-value was calculated as 2:95 · 10 25 at a = 0.05. In the most highly associated region, this analysis was repeated for windows of 20, 10 and 6 SNPs in sliding windows overlapping by n 2 1 SNPs in order to fine map the associated regions. This was carried out from approximately 5MB before and after the significant region.
Accounting for sample size difference between males and females: Sample sizes within this dataset are markedly different between males and females (see above and Table 1). A consequence of this may be that there is lower power to detect associations with male recombination rate. We repeated the heritability and GWAS analyses in sampled datasets of the same size within each sex. Briefly, 482 recombination rate measures (representing the total number in males) were sampled with replacement within the male and female datasets, and the animal model and GWAS analyses were repeated in the sampled dataset. This process was repeated 100 times, with sampling carried out in R v3.3.2. The observed and simulated heritabilities compared to see how often a similar results would be obtained. This was repeated for association at the most highly associated GWAS SNPs and regional heritability regions. The differences between the mean simulated values in each sex were investigated using a Welch two-sample t-test assuming unequal variances.
Haplotyping and effect size estimation: Haplotype construction was carried out to examine haplotype variation within regions significantly associated with recombination rate variation in the regional heritability analysis. SNP data from deer linkage group 12 was phased using SHAPEIT v2.r837 (Delaneau et al. 2013), specifying the linkage map positions and recombination rates for each locus. This analysis used pedigree information with the -duohmm flag to allow the use of pedigree information in the phasing process (O'Connell et al. 2014). Haplotypes were then extracted for the most significant window from the regional heritability analysis.
Effect sizes on ACC for the top GWAS SNPs were estimated using animal models in ASReml-R; SNP genotype was fit as a fixed factor, with pedigree relatedness fit as a random effect to account for the remaining additive genetic variance. To determine the effect sizes on ACC for the regional heritability analysis, animal models were run as follows: for a given haplotype, A, its effect was estimated relative to all other haplotypes combined, i.e., treating them as a single allele, B, by fitting genotypes A/A, A/B and B/B as a fixed factor. This was repeated for each haplotype allele where more than 10 copies were present in the full dataset.

Data availability
Raw data are publicly archived at doi: https://doi.org/10.6084/m9. figshare.5002562 . Code for the analysis is archived at https://github.com/susjoh/Deer_Recombination_GWAS. Supplementary data files are archived in the GSA figshare portal. Table S2 contains raw and quality controlled ACCs per individual/chromosome/meiosis. Table  S3 contains the ACC GWAS results for both sexes and in males and females only. Table S4 contains the ACC regional heritability results in both sexes and in males and females only. Table S5 contains the ACC regional heritability results in the most highly associated region on deer linkage group CEL12. Supplemental material available at Figshare: https:// doi.org/10.25387/g3.6166364.

RESULTS
Variation and heritability in autosomal crossover count Autosomal crossover count (ACC) was significantly higher in females than in males, where females had 4.32 6 0.41 more crossovers per gamete (Z = 10.57, P Wald ,0.001; Figure 1; Table 1); there was no effect of FID age or inbreeding on ACC (P . 0.05, Table S1). Females had significantly higher phenotypic variance in ACC than males (V P = 32.02 and 15.33, respectively; Table 1). ACC was heritable in both sexes combined (h 2 = 0.13, SE = 0.05, P = 0.002) and within females only (h 2 = 0.11, SE = 0.06, P = 0.033), but was not heritable in males alone (P . 0.05; Table 1). The remaining phenotypic variance was explained by the residual error term, and there was no variance explained by the n Table 1 Data set information and animal model results for autosomal crossover count (ACC). Numbers in parentheses are the standard error, except for Mean, which is the standard deviation. N OBS ; N FID and N xovers are the number of ACC measures, the number of focal individuals (FIDS) and the total number of crossovers in the dataset. The mean ACC was calculated from the raw data. V P and V A are the phenotypic variance and additive genetic variance, respectively. h 2 ; pe 2 and e 2 are the narrow-sense heritability, the permanent environment effect, and the residual effect, respectively; all are calculated as the proportion of V P that they explain. The additive genetic components were modeled using genomic relatedness matrices. Pðh 2 Þ is the significance of the V A term in the model as determined using a likelihood ratio test   Table S3 and sample sizes are given in Table 1.
n Table 2 The top five most significant hits from a genome-wide association study of ACC in (  correlation (r a ) between males and females was 0.346, but was not significantly different from zero or one (P LRT .0.05). This may be due to the relatively small sample size of this dataset resulting in a large standard error around the r A estimate, or the fact that ACC was not heritable in males. Sampling of 482 measures from each sex showed no difference in heritability estimates between the sexes, indicating reduced power to quantify heritable variation in the smaller male dataset (t = 0.242, P = 0.810, Figure S3). Raw and quality controlled ACC counts and positions for each FID and offspring combination per chromosome are provided in Table S2.

Genetic architecture of autosomal crossover count
Genome-wide association study: No SNPs were significantly associated with ACC at the genome-wide level (Figure 2, Table 2 and Table S3). The most highly associated SNP in both sexes was cela1_red_10_26005249 on deer linkage group 12 (CEL12), corresponding to position 26,005,249 on cattle chromosome 10 (BTA10). This marker was also the most highly associated SNP when considering recombination in trans, indicating that this region affects ACC across the genome (Table S3). The observed association was primarily driven by female ACC (Table 2, Figure 2). In females, the most highly associated SNP was cela1_red_10_25661750 on CEL12, corresponding to position 25,661,750 on BTA10. For both SNPs, sampling of 482 measures from each sex showed that the observed associations were significantly higher in females than in males when considering the same sample size (cela1_red_10_25661750: t = 18.60, P , 0.001; cela1_red_10_26005249: t = 4.89, P , 0.001; Figure S4). Based on its position relative to the cattle genome, cela1_red_10_26005249 was 600bp upstream of an olfactory receptor OR5AU1 and 24kb downstream from a gene of unknown function (ENSBTAG00000011396). There were four candidate genes within 1Mb of both loci, including TOX4, CHD8, SUPT16H and CCNB1IP1 (Figure 4; see Discussion). Similar results were obtained when considering recombination rate on all chromosomes excluding that on which the SNP occurred, indicating that all associations affect recombination rate variation in trans across the genome (Table S3).
Regional heritability analysis: The genome-wide regional heritability analysis of ACC showed a significant association in both sexes and in females only with a 2.94Mb region on CEL12 (Figure 3, Table 3). The most highly associated window (1.36 Mb) within this region contained 42 genes, including REC8 meiotic recombination protein (REC8; 20,810,610 -20,817,662 bp on BTA10). Detailed examination of this region in sliding windows of 6, 10 and 20 SNPs found the highest association at a 10 SNP window of 463kb containing 36 genes, including REC8 (Table 3). This region explained all heritable variation in recombination rate, with regional heritability estimates of 0.143  Table S4.
(SE = 0.053) and 0.146 (SE = 0.045) for all deer and females only, respectively. The sex-specific effect was supported by sampling of 482 measures, where females had consistently higher associations than in males (t = 19.03, P , 0.001, Figure S5). The total significant region after detailed examination was 3.01Mb wide, flanked by SNPs cela1_red_10_18871213 and cela1_red_10_21878407 ( Figure 4 & Table S5) and containing 87 genes. This wider region contained the protein coding region for ring finger protein 212B (RNF212B; 21,466,337 -21,494,696 bp on BTA10), a homolog of RNF212, which has been directly implicated in synapsis and crossing-over during meiosis in mice (Reynolds et al. 2013). Genetic variants at both RNF212B and RNF212 have been associated with recombination rate variation in humans, cattle and sheep (Kong et al. 2008;Ma et al. 2015;Johnston et al. 2016;Petit et al. 2017). While this region was close to the most highly associated SNPs from the genome-wide association study, there was no overlap between the two analyses, with the mostly highly associated regions separated by an estimated 5.5Mb (Figure 4). The mean r 2 LD between the top regional heritability window and the top GWAS SNPs was 0.258 for cela1_red_10_ 25661750 and 0.276 for cela1_red_10_26005249, with the top r 2 of 0.665 observed between the SNPs cela1_red_10_21807996 and cela1_ red_10_26005249 (Figure 4). A second region on linkage group 32 almost reached genome-wide significance in the regional heritability analysis, corresponding to the region 38.7 -41.3Mb on cattle chromosome 27. This region contained the locus topoisomerase (DNA) II beta (TOP2B); inhibitors of this gene lead to defects in chromosome segregation and heterochromatin condensation during meiosis I in mice, Drosophila melanogaster and Caenorhabditis elegans (Li et al. 2013;Gómez et al. 2014;Hughes and Hawley 2014;Jaramillo-Lambert et al. 2016; Figure 3, Table 3 and  Table S4). Full results for the regional heritability analyses are provided in Tables S4 and S5.
Effect size estimation: At the most highly associated GWAS SNP, cela1_red_10_26005249, carrying one or two copies of the G allele conferred 3.3 to 3.9 fewer crossovers per gamete in females (Wald P , 0.001) and 1.8 -2.8 fewer crossovers per gamete in males (P = 0.009; Table 4). The most highly associated SNP in females, cela1_red_ 10_25661750, had a significant effect on ACC in females (P , 0.001) but not in males (P . 0.05; Table 4). This locus conferred 2.03 more crossovers in A/G females and 13.68 more in G/G females; however, the latter category contained 7 unique measures in only two individuals, and so this estimate is likely to be subject to large sampling error.
A total of 17 haplotypes in the 10 SNP region spanning cela1_red_ 10_20476277 and cela1_red_10_20939342 had more than ten copies in unique individual females (Table S6). Of these, two haplotypes, AGGA-GAGAAG and AGAGAAGAGA, had a significant effect on ACC relative to all other haplotypes (Table 4 and Table S6, Figure S6). Haplotype AGGAGAGAAG increased female ACC by 2.4 crossovers per gamete in heterozygotes (P , 0.001); homozygotes for the haplotype were rare (13 measures in 4 individuals) and so the large effect size estimate was again likely to be subject to large sampling effects ( Table 4). The haplotype AGAGAAGAGA reduced female ACC by 2.2 crossovers per gamete in heterozygous individuals (P , 0.05; Table 4). The r 2 LD between haplotype AGGAGAGAAG and the two most highly associated GWAS SNPs was 0.464 and 0.885 for cela1_red_10_26005249 and cela1_red_10_25661750, respectively; for haplotype AGAGAAGAGA, it was 0.229 and 0.036 for cela1_red_10_26005249 and cela1_red_10_ 25661750, respectively.

DISCUSSION
In this study, we have shown that autosomal crossover count (ACC) is 1.2· higher in red deer females than in males, with females exhibiting higher phenotypic and additive genetic variance for this trait; ACC was not significantly heritable in males. Almost all genetic variation in females was explained by a 7Mb region on deer linkage group 12. This region contained several candidate genes, including RNF212B and REC8, which have previously been implicated in recombination rate variation in other mammal species, including humans, mice, cattle and sheep (Kong et al. 2008;Reynolds et al. 2013;Ma et al. 2015;Johnston n Table 3 The most significant hits from a regional heritability analysis of ACC in (A) Both sexes, (B) Females only and (C) Males only. Sliding windows were 20 SNPs wide with an overlap of 10 SNPs. Lines in italics are the most highly associated regions from detailed examination of significant regions -in each case these are for 10 SNP windows. The x 2 and P values are for likelihood ratio test comparisons between models with and without a genomic relatedness matrix for that window; values in bold type are significant the the genome-wide level. The SNP locus names indicate the position of the SNPs relative to the cattle genome assembly vBTA_vUMD_3.  Petit et al. 2017). Here, we discuss in detail the genetic architecture of individual recombination rate, candidate genes underlying heritable variation, sexual-dimorphism in this trait and its architecture, and the conclusions and implications of our findings for other studies of recombination in the wild.
The genetic architecture of individual recombination rate Using complementary trait mapping approaches, we identified a 7Mb region on deer linkage group 12 (homologous to cattle chromosome 10) associated with ACC. The most highly associated GWAS region occurred between 25.6 and 26Mb (relative to the cattle genome position), although this association was not significant at the genomewide level. The most highly associated regional heritability region occurred between 20.5 and 20.9MB, around 5Mb away from the top GWAS hits ( Figure 4); association at this region was significant at the genome-wide level and explained almost all of the heritable variation in ACC in both sexes and in females only. Most variation in mean ACC was attributed to two haplotypes within this region (Table 4 and Table S6; Figure S6).
At present, it is not clear why the results of the two analyses occur in close vicinity, yet do not overlap. Assuming homology with humans, cattle, sheep and mice (Ensembl release 91, Zerbino et al. 2018), the two regions are separated by the highly repetitive T-cell receptor alpha/delta variable (TRAV/DV) locus, which may contain up to 400 TRAV/DV genes in cattle (Reinink and Van Rhijn 2009;Figure 4). This region is of an unknown size in deer; relative to the cattle genome, these regions are separated by 4.72Mb, but the deer linkage map distance is estimated as 1.86 centiMorgans (cM). The sex-averaged genome-wide recombination rate in deer is 1.04cM/Mb, suggesting this genomic region may be shorter in deer  and that these two regions are in closer vicinity This is supported by both the linkage map distance and patterns of linkage disequilibrium between the associated loci, particularly at the associated haplotypes (see Results & Figure 4). In addition, the small sample size used in the current study may result in increased sensitivity to sampling effects and bias in the estimation of the relative contribution of SNPs to the trait mean (GWAS) or variance (Regional heritability). Further investigation with higher samples sizes, whole genome sequencing approaches and improved genome assembly may allow more accurate determination of the most likely candidate Figure 4 Detailed figure of genes, association statistics and linkage disequilibrium patterns at the most highly associated region on on CEL12 (homologous to BTA10) for all deer of both sexes. All X-axis positions are given relative to the cattle genome vBTA_vUMD_3.1. The top panel shows protein coding regions, with annotation for candidate loci. The central panel shows the results for the regional heritability analysis (where lines represents a sliding windows of 6, 10 and 20 SNPs with an overlap of n-1 SNPs) and the genome-wide association study (where points indicate single SNP associations). The dashed lines are the genome-wide significance thresholds (green = regional heritability, black = genome-wide association). The checked shaded area shows the position of the T cell receptor alpha/delta locus (see Discussion). Underlying data are provided in Tables S3 & S5. The lower panel shows linkage disequilibrium between each loci using allelic correlations (r 2 ). genes and potential causal mutations (coding or regulatory) within this species.
Candidate genes for recombination rate variation Regional heritability analysis: The most highly associated region in the regional heritability analysis contained the gene REC8, the protein of which is required for the separation of sister chromatids during meiosis (Parisi et al. 1999). It also contained RNF212B, a paralogue of RNF212. RNF212 has been associated with recombination rate variation in humans, cattle and sheep (Kong et al. 2008;Sandor et al. 2012;Johnston et al. 2016;Petit et al. 2017); the REC8/RNF212B region is associated with recombination rate in cattle, and has a large effect size on ACC phenotype than RNF212 in this species (Sandor et al. 2012;Ma et al. 2015). A second region on deer linkage group 32 almost reached genome-wide significance (Figure 3, Table 3 and Table S4). This region was relatively gene-poor, but contained 6 genes, including the candidate topoisomerase (DNA) II beta (TOP2B): inhibitors of this gene lead to defects in chromosome segregation and heterochromatin condensation during meiosis I in mice, Drosophila melanogaster and Caenorhabditis elegans (Li et al. 2013;Gómez et al. 2014;Hughes and Hawley, 2014;Jaramillo-Lambert et al. 2016). No association was observed at the region homologous to RNF212 (predicted to be at position 109.2Mb on cattle chromosome 6, corresponding to 57.576cM on deer linkage group 6) for the GWAS or regional heritability analysis.
Genome-wide association study (GWAS): Examination of annotated regions within 500kb of either side of the most significant GWAS SNPs identified three genes, TOX High Mobility Group Box Family Member 4 (TOX4), Chromodomain Helicase DNA Binding Protein 8 (CHD8) and SPT16 Homolog Facilitates Chromatin Remodelling Subunit (SUPT16H). These genes are involved in with chromatin binding and structure (SP16H, TOX4), histone binding (CHD8, SUPT16H), nucleosome organization (SP16H) and cell cycle transition (TOX4). One of these genes, SUPT16H, interacts with NIMA related kinase 9 (NEK9), which is involved with meiotic spindle organization, chromosome alignment and cell cycle progression in mice (Yang et al. 2012) and is a strong candidate locus for crossover interference in cattle (Wang et al. 2016). The SNP cela1_red_10_26005249 was 825kb from Cyclin B1 Interacting Protein 1 (CCNB1IP1), also known as Human Enhancer Of Invasion 10 (HEI10), which interacts with RNF212 to allow recombination to progress into crossing-over in mice (Qiao et al. 2014) and Arabidopsis (Chelysheva et al. 2012); this locus is also associated with recombination rate variation in humans (Kong et al. 2014).
Sexual dimorphism in genetic architecture of recombination rate: The results of this analysis suggest that there is sexual dimorphism in the genetic architecture of recombination rate variation in deer. Higher female recombination rates are typical in mammals (Brandvain and Coop, 2012) and in this system is driven by higher female recombination rates near centromeric regions . At present, the mechanisms driving differences in mean ACC within and between species are unknown (Lenormand and Dutheil, 2005;Stapley et al. 2017). Male ACC was not significantly heritable, although we could not rule out that this was a consequence of the smaller sample size relative to females ( Figure S3). No regions of the genome were significantly associated with male ACC in the regional heritability and GWAS analyses, but sampling did indicate that differences observed between male and female genomic associations were statistically significant ( Figures S4 & S5). Investigation of genetic correlations between males and females was inconclusive, as the r A of ACC was not significantly different from 0 or 1. The observed sex differences are consistent with previous studies of the genetic architecture of ACC in mammals, where a sexually-dimorphic architecture has been observed at the paralogous RNF212 region in humans and sheep (Kong et al. 2014; n Table 4 Effect sizes for the most highly associated GWAS SNPs and for the AGGAGAGAAG haplotype at the most highly associated regional heritability region. Models were run for each sex separately and included a pedigree relatedness as a random effect.  Johnston et al. 2016). Nevertheless, some observed associations were stronger when considering both male and female deer in the same analysis, for example at the most highly associated GWAS SNP, and the amplified signal for the regional heritability analysis on linkage group 33 (Figures 2 & 3), suggesting that there may be some degree of shared architecture within these regions.

Conclusions and implications for studies of recombination in the wild
We have shown that recombination rate is heritable in female red deer, and that it has a sexually dimorphic genetic architecture. Genomic regions associated with recombination rate in red deer are associated with this trait in other mammal species, supporting the idea that recombination rate variation has a conserved genetic architecture across distantly related taxa. A key motivation for this study is to compare how recombination rate and its genetic architecture is similar or different to that of model species that have experienced strong selection in their recent history, such as humans, cattle, mice and sheep. The heritability of recombination rate in deer was lower than that observed in other mammal systems (Dumont et al. 2009;Kong et al. 2014;Ma et al. 2015;Johnston et al. 2016), with no observed heritable variation present in male deer. While we were able to test their effects, we found no contribution of contribution of individual and common environmental effects on recombination rate (i.e., age, year of birth, year of gamete transmission); indeed, most phenotypic variance in recombination was attributed to residual effects. This suggests that despite some underlying genetic variation, recombination rate is mostly driven by stochastic effects, or otherwise unmeasured effects. This represents one of the smallest datasets in which recombination rate has been investigated, and so it may be that the observed effects are underestimated due to the small sample size, sampling effects, or perhaps that other genetic variants present in this species do not segregate in the Rùm deer population. Nevertheless, identification of clear candidate genes and their effects on phenotype represents a valuable contribution to understanding the genetic architecture of recombination more broadly. Ultimately, our findings allow future investigation of the fitness consequences of variation in recombination rate and the relationship between identified variants and individual life-history variation, to address questions on the maintenance of genetic variation for recombination rates, and the relative roles of selection, sexually antagonistic effects and stochastic processes in contemporary natural populations.