Mapping of quantitative trait loci for life history traits segregating within common frog populations

The evolution of complex traits is often shaped by adaptive divergence. However, very little is known about the number, effect size, and location of the genomic regions influencing the

variation of these traits in natural populations. Based on a dense linkage map of the common frog, Rana temporaria, we have localized, for the first time in amphibians, three significant and

nine suggestive quantitative trait loci (QTLs) for metabolic rate, growth rate, development time, and weight at metamorphosis, explaining 5.6–18.9% of the overall phenotypic variation in

each trait. We also found a potential pleiotropic QTL between development time and size at metamorphosis that, if confirmed, might underlie the previously reported genetic correlation

between these traits. Furthermore, we demonstrate that the genetic variation linked to fitness-related larval traits segregates within Rana temporaria populations. This study provides the

first insight into the genomic regions that affect larval life history traits in anurans, providing a valuable resource to delve further into the genomic basis of evolutionary change in

amphibians.

Revealing the genetic architecture behind adaptive processes is a fundamental issue in evolutionary biology. In particular, information about the location, effect size, and number of loci

controlling the life history-related traits is essential for understanding the mechanisms of evolutionary change (Lynch and Walsh 1998; Mackay 2001; Barton and Keightley 2002). In this

context, quantitative trait locus (QTL) mapping has been one of the most widely used tools to identify genomic regions that control important adaptive traits in wild populations (Slate 2005;

Charmantier et al. 2014; Bendesky et al. 2017).

Larval traits, such as developmental time, size at metamorphosis, and growth rate, are thought to be under strong natural selection (Collins 1975; Berven and Gill 1983; Houde 1997; Peckarsky

et al. 2001; Barton et al. 2014; Eck et al. 2015) driven by the environmental factors that affect growth opportunities such as predation, desiccation risk, and thermal conditions (Laurila

and Kujasalo 1999; Rodríguez-Muñoz et al. 2001; Fischer and Karl 2010; Tejedo et al. 2010). Among them, temperature is of paramount importance in shaping life-history variations (Atkinson

1994; Gotthard 2001; Angilletta et al. 2004). For instance, a low temperature slows down the metabolism, especially in ectotherms, with direct consequences on growth and development rates.

Ectotherms from high altitudes and latitudes often exhibit countergradient variations (i.e., rapid growth and development rates) in order to compensate for these unfavourable environmental

conditions and time constraints (reviewed by Conover et al. 2009). These organisms increase their food conversion efficiency and allocate their available energy at the expense of other

purposes (Angilletta et al. 2004); thus, a rapid growth might also entail some trade-offs related, for example, to locomotor performance (Cano and Nicieza 2006) or immune function (De Block

et al. 2008). Therefore, knowing the genomic basis of growth and metabolic rates is critical to understand the mechanisms behind the adaptive divergence in ectotherms. Furthermore, this is

highly relevant to understand how ectothermic organisms will deal with environmental variations such as climate change (Atkinson 1994; Angilletta et al. 2004; Umina et al. 2005; Johnston and

Bennett 2008).

Amphibians are ectothermic animals distributed worldwide and they inhabit a wide range of environmental conditions, which make them good models to study adaptive processes (Miaud and Merilä

2001; Beebee 2005). In addition, since planned crosses can be easily conducted by artificial fertilization, the quantitative genetic basis of larval life history traits can be characterized

(e.g., Berven and Gill 1983; Laurila et al. 2002; Palo et al. 2003; Laugen et al. 2005). These studies suggest that many larval traits have a significant heritable component. For instance,

development time, weight at metamorphosis, and GR exhibit significant additive genetic variation (Travis 1981; Berven 1987; Laurila et al. 2002; Palo et al. 2003; Knopp et al. 2007).

However, the magnitude of additive genetic variance and the strength and direction of their genetic correlations differ among populations and selective environments (Berven and Gill 1983;

Cano et al. 2004).

In addition to the increasing evidence about the heritable basis of life history traits, we also need complementary genetic information to uncover the molecular mechanisms responsible for

adaptive evolution (Conner and Hartl 2004). For instance, QTL studies inform about the genomic regions influencing phenotypic variations (Corva and Medrano 2001; Slate 2005; Beraldi et al.

2007; Rothschild et al. 2007; Lv et al. 2016).

Genomic regions controlling body weight and growth rate have been identified in a large number of endotherms such as mouse, pig, and chicken (e.g., Kerje et al. 2003; Jerez-Timaure et al.

2004) and, more recently, in ectotherms such as nematode and fish (e.g., Gutteling et al. 2007; Vasemägi et al. 2010; Lv et al. 2016). However, in amphibians, only development time has been

studied in relation to the evolution of paedomorphosis in Ambystoma (Voss et al. 2012; Page et al. 2013). The former study identified three QTLs (met1, met2, and met3) explaining 10–11% of

variation in development time by using a linkage map consisting of 185 molecular markers (Voss et al. 2012). The latter study showed that met1 genotype affected the expression of 200 genes

during larval development, linking this QTL with thyroid hormone signalling and mitochondrial energetics (Page et al. 2013). In other amphibians, no QTL mapping studies have been published

to date, although scans for footprints of selection (Bonin et al. 2006; Guo et al. 2016) and comparative transcriptome analyses (Yang et al. 2012) have identified a few candidate genes

potentially related to adaptation to high altitudes. In addition, the recent publication of a draft genome assembly and two dense linkage maps for the common frog, R. temporaria (Brelsford

et al. 2016; Palomar et al. 2017), allow to perform efficiently the genome-wide searches of QTLs and candidate genes.

Here, we used a high-density linkage map consisting of 7138 molecular markers to identify, for the first time, the genomic regions controlling the key larval life history traits in the

common frog, R. temporaria. We crossed two individuals from contrasting environments (i.e., high and low altitudes) and measured standard metabolic rate, growth rate, development time, and

weight at metamorphosis in their F1 offspring to estimate the effect size and number of QTLs associated to these traits.

A full-sib family was generated by crossing artificially the parents from two phenotypically and genotypically well-differentiated populations of R. temporaria (Choda 2014; see detailed

methods in Palomar et al. 2017). These populations are associated with different mitochondrial lineages and exposed to contrasting environmental conditions determined by elevation,

hydroperiod, and landscape structure. The male was captured from a mountain area, Vega de Candioches, León, (1687 m.a.s.l.) and the female from a lowland location, river Argonza valley, near

Bárcena Mayor, Cantabria, (551 m.a.s.l.), both in Northern Spain. Over 500 embryos were obtained by artificial fertilisation (eggs and sperm were obtained by pressing gently the animal

abdomens, a simple procedure during the breeding season) and maintained in dechlorinated water at 9 °C. At Gosner stage 25 (i.e., gill resorption completed and exogenous feeding started;

Gosner 1960), around 300 tadpoles were individualized in 0.8 l tanks with dechlorinated water following a fully randomised design. Larvae were fed with unrestricted rations of rabbit chow

(15% protein, 3% fat, 17% carbohydrate, 10% ash; Cargill España, Martorell, Barcelona, Spain) until metamorphosis. Animals were reared under conditions of constant photoperiod (12L:12D) and

temperature (14.0 ± 0.5 °C). Larvae were checked every day looking for metamorphs (Gosner stage 46). Metamorphs were euthanized with an overdose of Benzocaine (Ethyl 4-aminobenzoate; Sigma

Aldrich, ref.: 112909) and then frozen at −45 °C.

Traits measured in this study are connected with fitness and are related to each other. We measured standard metabolic rate, which is related to growth and development rates as well as to

the size and age at metamorphosis (Blackmer et al. 2005; Careau et al. 2008; Artacho and Nespolo 2009; Burton et al. 2011; Rosenfeld et al. 2015). Standard metabolic rate (SMR) was measured

at Gosner stage 33 (SD = 2.97). We used a flow-through respirometry system consisting of 24 cylindrical chambers (54 mm × 16 mm) immersed in water at a constant temperature of 14 ± 0.5 °C

(for a similar set up, see Álvarez and Nicieza 2005; Cano and Nicieza 2006). The metabolic chambers were supplied with oxygen-saturated water at a fixed flow rate of 700 ml/min controlled by

a 24-channel high-precision peristaltic pump (Model ISM934C; ISMATEC, Cole-Parmer GmbH, Germany). Larvae were unfed for 48h and then acclimated for 15h at the respirometry chambers prior to

measuring SMR. The tadpoles were kept in darkness and remained quiescent throughout SMR measurement. We measured the oxygen consumption in a flow-through system by using a thermostatted

cell (MC 100, Strathkelvin Instruments Ltd, Glasgow, UK) housing a microcathode oxygen electrode (Model 1302, Strathkelvin Instruments Ltd) connected to an oxygen meter (Model SI782

Single/Dual Channel Meter, Strathkelvin Instruments Ltd). The electrode was calibrated against air-saturated water (obtained from the header tank) and against a solution having zero oxygen

saturation (sodium sulphite in 0.01 sodium tetraborate). Bacterial oxygen consumption was prevented by using ultraviolet lamps. All the equipment was exposed to ultraviolet light for 30 min

before the procedure. For each tadpole, we measured the oxygen saturation at the outlet of a blank (empty) chamber and at the outlet of the tadpole chamber over a 5-min period. These

measurements were transferred via Strathkelvin software, 929 Oxygen System v01.02, and recorded in a computer for further analysis. SMR was calculated as follows:

where Vo2 (µg/h) is the rate of oxygen consumption, Vw is the flow rate (ml H2O/h) through the respirometry chamber, ΔCw is the difference in oxygen concentration between the blank and the

test chamber (µg O2/ml), and So2 is the solubility of oxygen in water (µg O2/ml) (Álvarez and Nicieza 2005; Cano and Nicieza 2006). Eventually, peaks in consumption derived from occasional

animal movements were identified and discarded. As a control for body-size variation, we used the residuals of the linear regression of SMR on tadpole mass.

In addition, we weighed the tadpoles weekly over a 4-week period (at 35, 42, 49, and 56 days after fertilization, relating to Gosner stages from 26 to 32) with a precision balance (±1 mg).

Since the increase of weight at these stages was linear (Appendix S1), the the growth rate (GR) was measured as the slope of the line that described the linear model between weight and time

for each tadpole. Individuals were also weighed at Gosner stages 42 (W42), emergence of the forelimbs, and 46 (W46), total reabsorption of the tail. Developmenttime (DT) was quantified as

the period between fertilization and Gosner stage 42.

DNA was extracted from 162 frozen metamorphs with the DNeasy Blood and Tissue Kit (Qiagen). Restriction site-associated DNA (RAD) libraries were prepared, genotypes were called, and a

linkage map was constructed as detailed in Palomar et al. (2017). Briefly, DNA from parents and progeny was digested with restriction enzymes PstI and BamHI and ligated to 94 modified

Illumina adapters with T4 DNA ligase. Agarose gel electrophoresis in an E-Gel® iBase™ Power System was used to do the size selection of 400-bp fragments. These fragments were amplified by

PCR, purified, and sequenced on two paired end lanes (2 × 100) with the Illumina HiSeq 2000 in an Illumina Genome Analyzer platform. After quality filter, 91-bp sequences were obtained. Sire

sequences were used to generate a reference to align the remaining individuals. PCR duplicates were deleted and sire sequences were clustered and assembled de novo to obtain this reference.

Parent sequences were aligned against this reference and filters of number of reads per contig (from 10 to 1000) and maximum mismatch (10%) were applied. RAD loci containing single

nucleotide polymorphisms (SNPs) fixed for alternative alleles in parents were discarded while loci containing heterozygous SNPs in each parent were used to align the progeny. Heterozygotes

were called when the minor allele count was >10%. In addition to SNPs, 113 microsatellite markers were amplified in a multiplex reaction as also detailed in Palomar et al. (2017). All

molecular markers segregating according to Mendelian fashion (X2 test, p value