Genome-wide variation in recombination rate in Eucalyptus

Variation between individuals

The difference in RR amongst parents (Fig. 1a) was highly significant (F_9,90 = 6.1, P < 0.001) and due mainly to the high RR of two F₁ parents (F1.4-F and F1.5-M) and low RR in the male and female parents of one of the F₂ families (F2.C). The RR of the F₁ families (3.1 ± 0.09 cM/Mb) was greater than that of the F₂ families (2.8 ± 0.11 cM/Mb), but not significantly so (F_1,8 = 3.9, P = 0.085). The F₂ families had higher levels of heterozygosity than the F₁ (2-tailed t ₈ = 7.04, P < 0.001). However while negative, the correlation between parental heterozygosity and RR was not significant. There was no significant differences in RR between sexes (F_1,8 = 0.39, P = 0.549). The genome-wide RR ranged between 2.71 and 3.51 cM/Mb between the 10 E. globulus parents, and averaged 2.98 cM/Mb (Table 1). These values were calculated across linkage maps ranging in size (i.e., genetic distance, estimated by Haldane’s mapping function) from 1442 to 1787 cM (Table 2), averaging 1545 cM.

Level and variable	Min	Max	Mean	CVe
Chromosome averages (n = 11)
Physical size (Mb)	39.02	80.10	55.10	0.29
Diversity (global HHw)b	0.196	0.229	0.196	0.04
Divergence (global F st )c	0.385	0.405	0.354	0.03
GC content (%)	35.60	38.29	37.23	0.02
RR (cM/Mb)	1.98	3.81	2.98	0.20
Individual genome-wide averages (n = 10)
Genetic distance (cM)	1442	1787	1545	0.08
Physical size (Mb)	508.3	569.5	543.1	0.03
Genome coverage (%)	84.13	93.55	89.37	0.03
Heterozygosityd	630	815	708	0.10
RR (cM/Mb)	2.71	3.51	2.98	0.08

^aDiversity and divergence were estimated by Hudson et al. [34]. The remaining variables were estimated from our E. globulus data and the E. grandis genome [19]

^bThe diversity within 6 Eucalyptus species, including E. globulus

^cDivergence between the 6 Eucalyptus species

^dEstimated from the total number of polymorphic markers derived from the DArT assay (see Methods)

^eCoefficient of variation

		Scaffold ID	1	2	3	4	5	6	7	8	9	10	11
Cross	Pa	Scaffold Size (Mb)b	40.3	64.2	80.1	42.0	74.7	53.9	52.4	74.3	39.0	39.3	45.5
F1.1	F	Interval Number	8	8	15	6	17	10	12	11	10	5	8
		Total length (cM)	122	170	174	99	125	199	147	206	132	91	122
		Phys. Size (Mb)	36.1	58.0	67.6	40.6	71.0	52.6	51.7	64.7	37.2	26.5	41.5
		Mean RR (cM/Mb)	3.38	2.93	2.57	2.44	1.76	3.78	2.84	3.19	3.55	3.44	2.94
	M	Interval Number	7	10	16	8	11	10	9	14	8	4	9
		Total length (cM)	150	176	159	82	129	173	105	128	99	113	139
		Phys. Size (Mb)	38.1	54.8	73.0	34.1	72.3	49.2	41.8	56.4	35.1	29.1	42.3
		Mean RR (cM/Mb)	3.93	3.21	2.18	2.40	1.78	3.52	2.51	2.27	2.82	3.88	3.29
F1.4	F	Interval Number	8	9	19	6	16	9	13	13	8	3	12
		Total length (cM)	176	199	184	98	145	183	180	224	118	117	163
		Phys. Size (Mb)	37.4	55.3	75.8	36.0	72.3	45.3	49.9	59.3	34.3	25.4	43.9
		Mean RR	4.71	3.60	2.43	2.72	2.01	4.04	3.61	3.77	3.44	4.61	3.72
	M	Interval Number	9	10	17	4	12	8	12	16	6	7	16
		Total length (cM)	153	144	142	110	143	154	115	177	82	99	137
		Phys. Size (Mb)	39.3	51.8	74.5	34.9	66.3	38.9	51.2	71.9	33.5	31.0	41.5
		Mean RR (cM/Mb)	3.90	2.78	1.91	3.15	2.16	3.96	2.25	2.46	2.45	3.19	3.30
F1.5	F	Interval Number	7	7	17	7	10	9	11	14	9	9	11
		Total length (cM)	115	142	136	101	109	165	152	170	104	104	144
		Phys. Size (Mb)	34.4	47.9	67.9	35.6	59.4	40.2	51.9	65.3	36.1	27.6	42.1
		Mean RR (cM/Mb)	3.35	2.96	2.00	2.84	1.84	4.11	2.93	2.60	2.88	3.77	3.42
	M	Interval Number	6	14	12	7	12	11	8	10	7	11	11
		Total length (cM)	139	157	184	116	166	181	173	206	109	136	151
		Phys. Size (Mb)	31.6	56.8	71.6	35.7	71.0	52.6	51.1	66.4	36.7	32.1	42.3
		Mean RR (cM/Mb)	4.40	2.76	2.57	3.25	2.34	3.44	3.38	3.10	2.97	4.24	3.57
F2.J	F	Interval Number	6	13	17	8	14	11	10	15	9	10	10
		Total length (cM)	107	167	139	101	129	163	101	158	105	162	168
		Phys. Size (Mb)	36.1	56.0	76.1	36.7	69.7	43.3	47.2	66.4	34.8	36.7	43.4
		Mean RR (cM/Mb)	2.96	2.98	1.83	2.75	1.85	3.76	2.14	2.38	3.02	4.42	3.87
	M	Interval Number	10	9	17	12	13	14	14	15	8	14	16
		Total length (cM)	129	165	153	92	141	164	142	181	121	130	124
		Phys. Size (Mb)	36.1	55.0	75.4	32.3	72.9	53.1	52.3	71.7	37.5	33.2	42.7
		Mean RR (cM/Mb)	3.57	3.00	2.03	2.85	1.93	3.09	2.71	2.52	3.23	3.92	2.90
F2.C	F	Interval Number	5	13	16	9	14	15	10	14	11	9	8
		Total length (cM)	105	164	160	110	144	168	133	137	117	120	118
		Phys. Size (Mb)	35.0	61.1	69.7	38.0	67.0	52.2	51.1	64.7	38.0	35.0	42.3
		Mean RR (cM/Mb)	3.00	2.68	2.30	2.90	2.15	3.22	2.60	2.12	3.08	3.43	2.79
	M	Interval Number	6	11	18	11	12	12	14	17	8	8	12
		Total length (cM)	119	172	139	92	136	189	121	178	102	113	124
		Phys. Size (Mb)	32.6	61.5	76.2	39.1	69.4	52.0	50.3	72.3	38.3	35.6	42.3
		Mean RR (cM/Mb)	3.65	2.80	1.82	2.35	1.96	3.64	2.41	2.46	2.66	3.18	2.93

^aP = Parent, F = Female, M = Male ^bPhys. Size = Physical Size

Variation between chromosomes

The 11 chromosomes of E. globulus differed markedly in RR. Averaged across parents, chromosome-wide RR ranged from 1.98 to 3.81 cM/Mb (Table 1). The difference between chromosome averages was highly significant when tested using parents as replicates (F_10,90 = 34.0, P < 0.001), and was also significant (F_10,89 = 92.7, P < 0.001) after removing the correlated effect of chromosome size (see discussion below). The variation between chromosomes was mainly due to the high RR in chromosomes 10, 1 and 6 and low RR in chromosomes 3 and 5 (Fig. 1b). The difference between chromosomes was moderately to highly positively correlated with previously published estimates for E. grandis x urophylla consensus maps ([26] n = 11, Pearson r = 0.90, P < 0.001; [25] n = 11, Pearson r = 0.92, P < 0.001) as well as parental maps [27] of E. grandis (n = 11, Pearson r = 0.72, P = 0.012) and E. urophylla (n = 11, Pearson r = 0.55) although the latter was not significant (P = 0.081). In no case was the interaction of chromosome with either sex (F_10,80 = 1.17, P = 0.327) or cross type (F_10,80 = 1.03, P = 0.426) significant, indicating that the chromosomal differences in RR were stable across sex and cross-type in this study.

Chromosomal variation in RR was also associated with differences in marker diversity (assuming Hard-Weinberg equilibrium [H _Hw ]) and species divergence (measured as the proportion of the genetic variance in a sub population relative to the total genetic variance [F _st ]). It was highly negatively correlated with the pooled genetic diversity (global H _Hw ) within the six Eucalyptus species studied by Hudson et al. [34] (n = 11, r = −0.75, P = 0.008). The correlation between chromosomal RR and H _Hw of each individual species, including E. globulus, were similar in size and all also negative (r = −0.60 to −0.74). The trends in global diversity were also evident regardless of whether markers were intra- or intergenic. However, for global H _Hw the correlation was only significant for intragenic markers (r = −0.65, P = 0.032). The global F _st amongst the six species studied by Hudson et al. [34] was positively correlated with RR at the chromosome level (r = 0.90, P < 0.001). This relationship with F _st was also found for both intragenic (r = 0.93, P < 0.001) and intergenic (r = 0.75 P = 0.008) markers.

Variation between individuals

The variation we found in genome-wide RR between individuals is consistent with reports in other studies [9, 37, 38] and may be due to environmental or genetic factors. RR has been shown to be affected by the environment in both plants and animals [39–41]. However, the role of the environment cannot be clarified in the present study as our parents come from different environments as well as being genetically distinct. Genetic factors which may influence RR, include sex and cross type. Sex specific differences have been reported in some taxa, with elevated RR often occurring in the gender subject to gametic selection [42]. However, no difference in RR between sexes was evident in E. globulus. There was a signal in our data, albeit not significant, that cross type could explain some of the variation between individuals, with lower average RR in the F₂’s compared with the F₁’s. However further study will be required to validate this. Irrespective of whether the genome-wide variation in RR we observed between individuals is due to genetic or environmental factors, it is clear that representative maps of the RR landscape in Eucalyptus species will need to integrate results from multiple crosses and individuals.

Our genome-wide average RR for E. globulus (2.98 cM/Mb) was higher than estimates from E. grandis (2.17 cM/Mb), E. urophylla (2.24 cM/Mb) [27], and E. grandis/urophylla consensus maps (1.58 cM/Mb and 1.95 cM/Mb; [25] and [26], respectively). However, differences in the scale over which RR was estimated is likely to have influenced the different estimates. Differences in genome coverage also no doubt influenced the estimates of RR between studies, particularly in the case of Kullan et al. [25], which estimated RR based on 152 windows of approximately 1 cM across the genome. Direct comparison of the estimates from different studies is also problematic because different mapping techniques produce different genetic distances and thus different estimates of genetic to physical distance. Our estimates of genetic distance (cM) were derived from Joinmap software, using the maximum likelihood algorithm, which uses Haldane’s mapping function. The estimates of genetic distance used by Petroli et al. [26] were derived from RECORD software [43]. RECORD produces distances very similar to the Joinmap maximum likelihood algorithm, if Haldane’s mapping function is used [44]. However, Silva-junior and Grattapaglia [27] used Kosambi’s mapping function, which produces smaller estimates of genetic distance relative to Haldane’s function. Similarly, both Kullan et al. [25] and Petroli et al. [26] used the regression algorithm and Kosambi’s mapping function in Joinmap, which produce considerably smaller genetic distances than Joinmap’s maximum likelihood algorithm. Therefore, the greater estimate of genome-wide RR in our study may simply reflect greater estimates of genetic distance due to methodological issues. Clearly, further study is required to draw robust conclusions regarding genome-wide differences in RR within and between eucalypt species, ideally using linkage maps constructed from similar techniques in multiple representatives of each species.

Variation between chromosomes

There was substantial variation in RR between chromosomes in our study, indeed chromosome was the most significant explanatory variable in our analysis of variance. Few studies in plants have examined genome-wide RR using multiple crosses, which are needed to provide the replication to test for significant differences between chromosomes. Nonetheless, chromosomal variation in RR has been demonstrated in a range of taxa [45, 46], including plants such as maize in which intraspecific variation was of a similar magnitude to our findings [14]. Chromosomal variation in RR is likely to at least in part reflect variation in chromosome size, due to the requirement for at least one crossover per chromosome per meiosis [45]. However, in our study the variation in RR between chromosomes remained significant after removing the correlated effect of chromosome size, arguing that other aspects of genomic architecture also contribute to this variation.

Despite the significant differences in genome-wide RR between individuals, three lines of evidence suggest that the relative RR in each chromosome is highly stable; not only within E. globulus, but between Eucalyptus species, and across evolutionary time. First, the differences between chromosomes in E. globulus were evident even when the significant effect of variation amongst the individuals was removed; and there was no significant interaction between chromosome and either sex or cross type within E. globulus. Indeed, the difference in RR amongst chromosomes was nearly two fold, whereas that amongst individuals was only 1.3 fold. Second, the high positive correlation between our chromosomal estimates of RR from E. globulus (section Maidenaria) and those in E. grandis/urophylla (section Latoangulatae) [25, 26] suggests the relative chromosomal differences in RR are stable across these Eucalyptus sections, despite an estimated 10 to15 Ma since their divergence [47]. Third, a positive correlation between GC content and RR is expected based on studies in other taxa (see below). The high positive chromosome-level correlation between GC content (from the E. grandis genome) and RR (measured from E. globulus map distance against the E. grandis genome), also points to stability of RR between E. grandis and E. globulus. Further, because our estimate of RR is from a single generation of recombination, while GC content in E. grandis reflects historical RR, the high positive correlation also points to a long term stability of RR at the chromosome level.

Association with genomic features and population genetic parameters

The chromosome level correlation between GC content and RR we observe is consistent with the theory of ‘GC-biased gene conversion’, whereby DNA repair during recombination promotes enrichment in GC bases, accounting for the positive correlation between RR and GC content seen in most eukaryotes studied to date [9, 48]. The particularly strong positive correlation observed between chromosome level RR and GC in E. globulus may reflect the unusually high GC content seen in the family Myrtales, compared with other eudicots [48]. In contrast, Arabidopsis is relatively GC poor and GC content and RR are uncorrelated, in line with a general pattern of greater GC enrichment in more derived plant genomes [48].

Our results also suggest a strong positive correlation between chromosomal estimates of RR and gene density. This finding is consistent with the positive correlation between RR and gene density observed at the whole-genome level in a range of taxa [9]. In plants such as maize, rice, wheat, and A. thaliana much of this correlation is attributable to the low density of genes in recombination-poor heterochromatin [49]. While in vertebrate genomes, Nam and Ellegren [7] hypothesise neutral processes associated with recombination, including a bias in deletions over insertions, contribute to a more compact genome structure and greater gene density in regions of elevated RR. However, the correlation between gene density and RR varies in strength and direction amongst taxa [9, 36, 50]. The positive correlation we observe at the chromosome level is particularly strong relative to other reports. Consistent with our findings, gene density was positively correlated with historical recombination rate in E. grandis (i.e., population scaled recombination rate, measured in 100 kb windows across the genome; [27]), suggesting the relationship between gene density and RR may be widespread in Eucalyptus, at least across the subgenus Symphyomyrtus.

A negative correlation between RR and some categories of transposable elements was also evident in our study. Transposable elements comprise over 50 % of most eukaryotes genomes and have a significant impact on genome evolution [51, 52]. In particular, it is well accepted that the accumulation of transposable elements is a major factor driving genome expansion [53]. The relative density of transposable elements across the genome reflects a balance between transposable element insertion and deletion. Transposable element deletion is believed to be primarily driven by recombinational processes, consistent with the observed transposable element accumulation in areas of low recombination in many taxa [49, 54] and in line with the correlation we observe. The mechanisms contributing to transposable elements accumulation in regions of low RR are the subject of some debate. Leading hypothesis involve relaxation of various forms of selection in these regions (see Dolgin and Charlesworth [55]); biased insertion of transposable elements in areas of low RR [56]; and a general deletion bias associated with recombination [7]. Regardless of the causative factors, in our study it is possible that in chromosomes with relatively low RR, a greater density of transposable elements has contributed to a reduced gene density, and thus contributed to the positive correlation we observe between RR and gene density at the chromosome level.

A notable finding to emerge from this study is our strong negative correlation between RR at the chromosome-level and species-wide genetic diversity. The common trend in animals is for a positive association between RR and genetic diversity, whereas the association in plants is more variable, but generally zero or positive [36]. The relationship between RR and genetic diversity may be affected by the scale and type of measures of RR (e.g., chromosome averages versus smaller genomic regions) and genetic diversity (i.e., within individuals, populations or species) both of which differ widely between studies [9]. Nevertheless, only in Oryza species (wild and domesticated rice) have negative correlations between RR and genetic diversity been previously reported [36]. In O. sativa, this association involved 100 kb intervals across three chromosomes (r = −0.183 and −0.314; for O. sativa ssp. indica and ssp. japonica, respectively) [57].

The more commonly observed positive correlation between RR and genetic diversity is often attributed to the impact of linked selection [9, 36]. Linked selection predicts a positive relationship between RR and genetic diversity due to the homogenising effects of selection being spread further in regions of low recombination [58]. However, as noted above there is disparity in this relationship between species [36] suggesting other factors also influence how selection impacts the genomic distribution of diversity, such as gene density [59]. Gene dense regions will be on average subject to greater selection, reducing diversity. Therefore, a more general prediction for the effect of selection across genomes is that the distribution of genetic diversity will be related to the density of elements under selection, i.e., gene density (assuming genes are the most likely targets of selection) [59]. This was the case in a recent study of E. grandis, where nucleotide diversity (at the population level) was negatively correlated with gene density in 100 kb windows across the E. grandis genome; and amongst the genomic features analysed, gene density yielded the highest correlation with nucleotide diversity, far exceeding the strength of the correlation between historical recombination rate and nucleotide diversity [27].

Thus when gene density and RR are strongly positively correlated, as in E. globulus, this may counter the influence of RR on linked selection, disrupting the relationship between RR and genetic diversity [36]. Indeed, there is growing evidence that the effects of linked selection are not confined to areas of low recombination but are often pervasive throughout genomes [36, 60], including in the forest tree poplar [61]. In support of this hypothesis, many plants exhibit a positive correlation between RR and gene density and a weak to no correlation between RR and genetic diversity [36, 60, 62]. In the case of Oryza sativa [57] and E. globulus the strong negative correlation between RR and diversity is likely to reflect, at least in part, the strong correlation between gene density and RR. However, in the case of eucalypts the contribution of other factors, such as their high genetic diversity [19, 27, 29] (see introduction) other facets of genome architecture, and the history of population size within each species, requires further investigation. Further, we cannot completely dismiss the possibly that other mechanisms are involved. For example, diversity per se could be affecting RR. It is well known that sequence divergence between individuals can inhibit recombination in adjacent regions of the genome, associated with the DNA mis-match repair system [63, 64] and this process could also contribute to the negative correlation between RR and diversity we observe.

Linkage map construction and comparison with physical map

Recombination rate (RR) was estimated in 10 unrelated genotypes, which were the parents of 5 full-sib pedigrees. These were all inter-provenance crosses and included three F₁ (each with 183 to 184 individuals) and two ‘outbred F₂’ (172 and 503 individuals) pedigrees. Together, these pedigrees sampled a diverse section of the natural distribution of E. globulus (see Hudson et al. [33], Freeman et al. [65]). All of the pedigrees had been previously used for linkage map construction [33, 65] and the progenies genotyped with diversity array technology ([DArT], DArT P/L Ltd Canberra, Australia) [66, 67] and microsatellite markers as described in Hudson et al. [33] and Freeman et al. [65]. To estimate RR in this study, we firstly rebuilt each of the 10 parental linkage maps using a subset of the markers used in the original studies. This was done to maximise the accuracy of RR estimates from the genome-anchored markers. The subset of markers were selected based on having a known and unique physical location on the BRASUZ1 reference genome of E. grandis [19] based on basic local alignment search tool (BLAST) searches, and linkage mapping criteria as described below.

For BLAST searches, adaptors were removed from the DArT sequences using Cutadapt V1.2.1 [68], then BlastN searches were performed in the Blastall V2.2.26 suite (alignment length > = 50 bp; e-value = 10⁻⁵⁰). Two BLAST analyses were conducted. Firstly, an auto-blast amongst all marker sequences was performed to identify redundant markers. Secondly, the physical location of non-redundant markers was determined through a BLAST search against the reference genome. For the latter, BLAST results were classified according to the criteria developed by Salse et al. [69]: the cumulative percentage of sequence identity (CIP) and the cumulative alignment percentage (CALP), based on the length of high-scoring segment pairs (HSP) and the sum of all HSP lengths (AL). For each marker, the most plausible position on the BRASUZ1 reference genome was identified based on the CIP value with a maximum CALP (markers with CALP > 200 % were removed) which minimised RR.

Of the non-redundant markers anchored to unique genome positions, four main mapping criteria were used to select the markers used to estimate RR. Ordered by priority, these criteria were: i) genome coverage; ii) segregation; iii) linkage statistics; and iv) RR. Markers with weak linkage statistics were excluded, with the exception of a few markers which were included to maximise genome coverage, provided their map order matched the physical map. Uniparentally inherited markers segregating 1:1 were selected preferentially to biparentally inherited markers segregating 3:1 (preventing auto-correlations when comparing recombination estimates in parental maps from the same pedigree), although markers segregating 3:1 were selected in cases where they increased map coverage. Linkage analysis was performed using the maximum likelihood algorithm in Join-Map v.4 [70]. Map distance in cM was calculated using Haldane’s mapping function. To estimate RR (cM/Mb), we aligned linkage groups (LGs) against the 11 main scaffolds of BRASUZ1. MapChart 2.2 software [71] was used to compare physical and genetic maps. Markers were accepted when the RR between adjacent markers was less than 25 cM/Mb, since greater RR probably reflects genotyping error(s).

Variation in RR and relationships with genome features and population genetic parameters

RR was estimated by comparing 2,452 meioses (i.e., 1,226 progenies in total across the 5 pedigrees × 2 parents in each pedigree; see above) in E. globulus over the 640 Mb of the reference genome of E. grandis. Variation in RR was studied at two different levels: within chromosomes and genome-wide. At the chromosome level, two estimates of RR were used: i) a RR for each chromosome of each parent, calculated as the ratio of the total LG length in cM to the corresponding physical size per scaffold in Mb; and ii) an average chromosome RR (n = 11) across all 10 parents. The RR for each chromosome of each parent (n = 110; Table 2) was used to test the significance of the parent and chromosome main effects, using a fixed effect model with no interaction and Type 3 Sum of Squares. The chromosome averages (n =11) were used for studying the (Pearson’s) correlation with other chromosomes features (chromosome average data for the main features are presented in Table 4). These analyses were undertaken using PROC MIXED and PROC CORR of SAS respectively. At the genome-wide level, for each parent, we used the total genetic distance and the corresponding physical size, to estimate an average RR.

Heterozygosity, diversity and divergence indexes were correlated against RR at the chromosome level. Relative heterozygosity was estimated for each individual and at the chromosome level for each individual, from the total number of polymorphic markers derived from the DArT assay. Diversity and divergence indexes were derived from independent samples [34]. These diversity indexes were calculated from six Eucalyptus species from subgenus Symphyomyrtus, using ~ 90 range-wide samples per species. The six species belong to three sections within Symphyomyrtus: section Latoangulatae (E. grandis and E. urophylla), section Maidenaria (E. globulus, E. nitens and E. dunnii) and section Exsertaria (E. camaldulensis). These species are among the most important hardwood plantations species world-wide and are representative of the sections (lineages) that contain the other economically important eucalypt species [29]. We used the average diversity (H _Hw ) in each species, as well as the average diversity within, and divergence (F _st ) between, the six species in our correlation analyses at the chromosome level. The ‘global’ estimates (i.e., those averaged across six species) of diversity and divergence were also analysed after partitioning the markers into those within genes and those >5 kb from a gene (i.e., intra- versus intergenic markers; following [34]).

We also tested for relationships between RR in E. globulus and genomic features of the Eucalyptus reference genome [19], based on our estimates of average RR for each chromosome, using Pearson’s correlation. These chromosomal features included: physical size of each E. grandis chromosome scaffold (in Mb), the GC content, average gene density, the number of tandem duplicated genes and the density of transposable elements (also broken into subcategories of DNA transposable elements, retrotransposons and uncategorised elements, and superfamilies within the DNA transposable elements and retrotransposon subcategories).