Animals and management
The experimental animals were Romane ewes reared at the INRAE La Fage Experimental Farm (Causse du Larzac, Saint-Jean Saint-Paul, southern France) between 2006 and 2019 (n = 1034) [35]. Ewes were reared exclusively outdoors on approximately 280 ha of rangelands, in a flock of 250 reproductive females present each year. The main management features of this farming system have been previously described in detail [36,37,38]. Briefly, the farming system was based on a productive flock reared exclusively in extensive harsh conditions while limiting supplementation, in order to investigate the capacity of ewes to fend for themselves. In the autumn, before mating began, dry ewes successively grazed native and fertilized rangelands (6% of the total surface). A single mating period took place at the end of the autumn and first mating occurred at 8 or 20 months of age depending on the ewe’s live weight and experiment. During the winter and the second half of pregnancy, ewes were gradually supplemented with conserved feedstuffs (i.e., hay and silage produced on the farm) and barley due to the absence of grazeable biomass on the rangelands (for the detailed composition of the diet, see [37]). Lambing took place outdoors in the spring (April) and ewes suckled lambs for approximately 80 ± 4 days while they successively grazed fertilized and native rangelands. During the summer, dry ewes grazed the senescent vegetation due to drought on large paddocks containing a high proportion of shrubs (up to 30%). The Romane ewes produced an average of 2.2 live lambs per lambing in our conditions.
Measurements
All the ewes were individually monitored for their body reserves (BR) through body condition score (BCS), body weight (BW), subcutaneous back fat (BF) and muscle (BM) depth, and their pedigree information was recorded. All the data were recorded in an INRAE experimental database for sheep and goats (GEEDOC). BCS measurements were performed according the original scale described by Russel et al. [9] (i.e., ranging from 1, emaciated, to 5, obese) and the scale was adapted with a subdivision of 0.1 increments instead of 0.25. The same two operators systematically recorded the BCS measurements over the 14-year period and underwent regular training sessions for calibration. Subcutaneous back fat and muscle depths were measured by a real-time ultrasound system on the 12th rib. The two technicians mentioned above were also involved in BF and BM measurements (performed over the last 4 years of the period) and underwent regular training. The measurements used in this study were collected on a regular basis during one to three productive cycles, according to the following physiological stage schedule: mating (M, 15 days before mating), early pregnancy (Pa, 39 ± 11 days after mating), two-thirds of pregnancy (Pb, 101 ± 11 days after mating), lambing (L), early suckling (Sa, 17 ± 10 days after lambing), middle of the suckling period (Sb, 42 ± 10 days after lambing), weaning (−W, 80 ± 10 days after lambing) and post-weaning period (Wp, 149 ± 11 days after lambing). To characterize body reserves dynamics (BRD), differences in BCS between pairs of physiological stages were calculated and analyzed (BCS-Pa:L, BCS-L:Sa, BCS-Pa:W, BCS-M:Pa, BCS-W:Wp and BCS-W:M; described by Macé et al. [24]). In addition to BCS, BW was also measured at all the stages described above, and BF and BM were measured at mating, lambing and weaning.
Descriptive statistics
Analyses of variance were carried out, taking the repeated measurements into account, using the MIXED procedure of the Statistical Analysis System (SAS version 9.4; SAS Institute Inc., Cary, NC, USA) to test relevant effects and interactions affecting phenotypes. The age at first lambing, the parity of the ewe, the litter size class and the year of measurement were identified as fixed effects and a random animal effect was included to consider repeated measurements on the same animal. The age at first lambing effect took account of ewes that lambed for the first time at 1 or 2 years of age (classes 1 and 2, respectively). The parity effect took account of first, second, third and more lambing (classes 1, 2 and 3, respectively). The effect of litter size class took account of the number of lambs born and suckled (class 0, empty ewes or lambing but without suckling lambs; class 1, ewes lambing and suckling singleton from L to W; class 2 ewes lambing more than singleton and suckling one lamb; class 3, ewes lambing and suckling twins; and class 4, ewes lambing and suckling more than twins). At mating, the effect of litter size class from the previous parity was considered. The first-order interactions between age at first lambing × litter size class and parity × litter size class were tested. An effect was considered significant if P < 0.05.
Phenotypes
First, raw phenotypes were adjusted for significant fixed effects consisting of age at first lambing, parity of the ewe, litter size class, and year of measurement fitting the following linear model:
$$\mathrm{y}=\mathrm{Xb}+\mathrm{e},\kern0.5em \left[\mathrm{Model}\ 1\right]$$
where y is the vector of observations for one of the BR traits; b is the vector of fixed effects; and e is the vector of random residuals with incidence matrix X. The residuals resulting for this model (1) were used as adjusted phenotypes to fit an animal model to estimate individual values using ASREML 3.0 software [39] with the following linear mixed model:
$$\begin{array}{cc}\mathrm y^\ast=\mathrm{Zc}+\mathrm e,&\left[\mathrm{Model}\;2\right]\end{array}$$
where y* is the vector of adjusted phenotypes; c is the vector of random animal effects; and e is the vector of random residuals with incidence matrix Z. c and e were assumed to be normally distributed with means equal to zero and (co) variances Iσ2c, Iσ2e, respectively. I are identity matrices of appropriate size. Model 2 was fitted with an identity matrix I for random animal effects fitted instead with the pedigree-based relationship matrix, as performed in Mucha et al. [40]. The resulting estimated animal values were used as phenotypes for subsequent genomic analyses instead of estimating breeding values. This was done to take repeated measurements for each individual into account, considering animal effect as a permanent environmental effect and not to use estimated breeding values that give very high false positive rates (i.e. increase type 1 error rates) due to influence of pedigree relationships (contributions of parents and relatives) on an animal’s value [41].
SNP genotypes and quality control
Sheep were commercially genotyped with Illumina Ovine SNP15K (i.e., 16,560 Single Nucleotide Polymorphisms (SNP); low density (LD)), SNP50K (i.e., 54,241 SNPs, medium density (MD)) or SNP600K (i.e., 606,006 SNPs; high density (HD)) beadchips. Genotypes were established as part of the research projects: “SheepSNPQTL”, “COMPAGNE”, “RomaneIteDomum”, “iSAGE” and “SMARTER”. A total of 1034 phenotyped female animals were genotyped and distributed as follow: 820 ewes were genotyped with the MD beadchip, 167 ewes were genotyped with the LD beadchip, and 47 ewes were genotyped with the HD beadchip. Among the phenotyped and genotyped females, 554 have their dam genotyped (i.e., 389 genotyped dams out of a total of 700 dams) and 965 have their sire genotyped (i.e., 49 genotyped sires out of a total of 60 sires). Dams and sires were genotyped with either MD or HD chips (374 parents and 64 parents, respectively). Concerning ewes genotyped with LD chip, both parents were genotyped.
Individuals with a call rate below 0.95 and with Mendelian inconsistencies were discarded (i.e., five genotyped and phenotyped ewes were removed). The SNPs were removed from further analyses if they were not in Hardy-Weinberg equilibrium, had a minor allele frequency below 1% or had a call rate below 0.98. PLINK software was used to detect incompatible genotypes between sires, dams and offspring [42]. When the total number of incompatible SNPs was more than 2% of all SNPs, ewes were kept in the analyses and the parents in error were replaced by missing values (14 ewes concerned). For ewes genotyped with the HD chip, only SNPs present on the MD chip and LD chip were kept. The ewes genotyped with the LD chip were imputed to the MD chip using Fimpute software, as were LD chip SNPs absent in the MD chip [43]. This resulted in a data set, used for QTL analyses, containing 1034 ewes with BR phenotypes genotyped for 48,593 autosomal SNPs.
QTL detection method
A genome-wide association study (GWAS) was performed using a univariate linear mixed model (LMM) to account for relatedness and population structure, as implemented in GEMMA v0.94.1 software [44], assessing significance with the Wald test. The statistical model used to test one marker at a time was:
$${\displaystyle \begin{array}{cc}{\mathrm{y}}^{\#}=\mathrm{Wp}+\mathrm{Zu}+\varepsilon, & \left[\mathrm{Model}\ 3\right]\end{array}}$$
where y# is the vector of phenotypes (i.e. solutions from vector c in model 2) for all individuals adjusted for fixed effects and accounted for repeated measures of animals; W is the vector of genotypes at the tested marker; p is the effect of the tested marker; u is a vector of random additive genetic effects distributed according to N(0, Aλτ− 1), where λ is the ratio of the additive genetic variance and the residual variance τ− 1, A is the additive relationship matrix, and Z is the incidence matrix (identity matrix in this case); ε is a vector of residuals distributed according to N(0, Iτ− 1), where I is the identity matrix.
GEMMA software implements the Genome-wide Efficient Mixed Model Association algorithm. The first step of the analyses included the estimation of the relatedness matrix. The resulting metrics were included in the second step (GWAS), allowing for the adjustment for both relatedness and population structure. The significance thresholds at the chromosome-wide level (BONFchr i) and genome-wide level (BONFgeno) were obtained using the Bonferroni method that accounts for multiple testing assuming that the number of independent tests was equal to the number of SNPs analyzed [45]. The formulas to obtain thresholds were the following:
$${\mathrm{BONF}}_{\mathrm{chr}\ \mathrm{i}}=-{\mathit{\log}}_{10}\left(\frac{\alpha }{SNPs\ i}\right)$$
$${\mathrm{BONF}}_{\mathrm{geno}}=-{\mathit{\log}}_{10}\left(\frac{\alpha }{SNPs\ geno}\right)$$
where SNPs i is the number of SNPs for chromosome i, and SNPs geno is the total number of SNPs at the genome level (i.e., 48,593 autosomal SNPs) and considering α = 5%. The resulting genome-wide significance threshold was equal to 5.94. GWAS was performed for all of the 26 autosomal chromosomes. The percentage of variance explained by each SNP was calculated as follows:
$$\%{\sigma}_p^2=100\left(\frac{2p\left(1-p\right){\alpha}^2}{\sigma_p^2}\right)$$
where \({\sigma}_p^2\) is the phenotypic variance of the trait (i.e. total phenotypic variance was obtained from variance components in model 2), p is the frequency of the allelic substitution effect of the SNP, and α is the estimated allelic substitution effect of the SNP [46]. When a given trait was significantly affected by multiple variants, the reported top SNPs are SNPs that have the highest –log 10 (P-value) among the significant SNPs in a 1-Mbp window.
The annotated candidate genes that were closest to the top SNPs were identified using the Ensembl release 104 of the sheep reference genome OAR v3.1 [47]. Position of the top SNP of QTLs and the annotated genes closest to the top SNPs were updated with the reference genome Rambouillet v1.0 (Additional file 2: Table S3).
LEPR structure and effect of LEPR genotype
Multiple LEPR protein sequence alignment of several mammals was performed using Weblogo software [23, 48]. Protein sequences are available at NCBI (Mus musculus NP_666258.2, Ratus norvegicus NP_036728.1, Homo Sapiens NP_002294.2, Sus scrofa NP_001019758.1, Bos taurus NP_001012285.2, Ovis aries NP_001009763.1 incomplete sequence) [49]. The complete LEPR protein sequence in sheep (W5PL31, 1165 amino acids) was obtained from the UniProt data base [50] from transcript ENSOART00000011314.1 (Ensembl), and we introduced the proline to serine substitution at position 1019.
Analyses of variance were performed to test the effect of the LEPR genotype (i.e., SNP oar3_OAR1_40,857,869) on BR phenotypes using the MIXED procedure of SAS (version 9.4; SAS Institute Inc., Cary, NC, USA). The dependent variables were phenotypes used in GWAS for BR and BRD at each physiological stage, as described above. The three possible LEPR genotypes were fitted as a fixed explanatory variable and a significance threshold of p < 0.05 was selected. The Varcomp procedure of SAS was used to fit the genotype effect as random and to estimate the proportion of variance explained by the genotype. The animal was included as a random effect in both the Mixed and Varcomp models.
The effects of the LEPR genotype on BW, BF and BM were independently analyzed at each physiological stage with a linear mixed model using the mixed procedure of SAS. The LEPR genotype was included as a fixed effect, whereas the animal was treated as a random effect. Appropriate following fixed effects, parity, litter size class, age at first lambing and year of measurement, were also including in mixed models depending on the trait analyzed.
