Front. Genet. Frontiers in Genetics Front. Genet. 1664-8021 Frontiers Media S.A. 745284 10.3389/fgene.2021.745284 Genetics Original Research Genetic Variations and Differential DNA Methylation to Face Contrasted Climates in Small Ruminants: An Analysis on Traditionally-Managed Sheep and Goats Denoyelle et al. Adaptation and Acclimation in Ruminants Denoyelle Laure 1 2 de Villemereuil Pierre 3 Boyer Frédéric 1 Khelifi Meidhi 1 Gaffet Clément 1 Alberto Florian 1 Benjelloun Badr 1 4 Pompanon François 1 * Univ. Grenoble Alpes, Univ. Savoie Mont Blanc, CNRS, LECA, Grenoble, France GenPhySE, Université de Toulouse, INRAE, ENVT, Castanet Tolosan, France Institut de Systématique, Évolution, Biodiversité (ISYEB), École Pratique des Hautes Études | PSL, MNHN, CNRS, SU, UA, Paris, France Institut National de la Recherche Agronomique Maroc (INRA-Maroc), Centre Régional de Beni Mellal, Beni Mellal, Morocco

Edited by: Xingbo Zhao, China Agricultural University, China

Reviewed by: Yu Jiang, Northwest A and F University, China

Juliane Friedrich, University of Edinburgh, United Kingdom

 *Correspondence: François Pompanon, francois.pompanon@univ-grenoble-alpes.fr

This article was submitted to Evolutionary and Population Genetics, a section of the journal Frontiers in Genetics

 28 09 2021 2021 12 745284 21 07 2021 02 09 2021 Copyright © 2021 Denoyelle, de Villemereuil, Boyer, Khelifi, Gaffet, Alberto, Benjelloun and Pompanon. 2021 Denoyelle, de Villemereuil, Boyer, Khelifi, Gaffet, Alberto, Benjelloun and Pompanon

This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

The way in which living organisms mobilize a combination of long-term adaptive mechanisms and short-term phenotypic plasticity to face environmental variations is still largely unknown. In the context of climate change, understanding the genetic and epigenetic bases for adaptation and plasticity is a major stake for preserving genomic resources and the resilience capacity of livestock populations. We characterized both epigenetic and genetic variations by contrasting 22 sheep and 21 goats from both sides of a climate gradient, focusing on free-ranging populations from Morocco. We produced for each individual Whole-Genome Sequence at 12X coverage and MeDIP-Seq data, to identify regions under selection and those differentially methylated. For both species, the analysis of genetic differences (FST) along the genome between animals from localities with high vs. low temperature annual variations detected candidate genes under selection in relation to environmental perception (5 genes), immunity (4 genes), reproduction (8 genes) and production (11 genes). Moreover, we found for each species one differentially methylated gene, namely AGPTA4 in goat and SLIT3 in sheep, which were both related, among other functions, to milk production and muscle development. In both sheep and goats, the comparison between genomic regions impacted by genetic and epigenetic variations suggests that climatic variations impacted similar biological pathways but different genes.

 adaptation methylated DNA whole genome Morocco small ruminant acclimation Introduction

The evolution and distribution of species is driven by the variation of their environment. Optimization of the individual’s phenotype to fit the environment may occur on the long-term, as populations adapt to the local values of environmental drivers (i.e., native individuals having on average a higher fitness than migrants, Savolainen et al., 2013). The access to whole genome sequences now allows to study the genomic bases of local adaptation by identifying genes and genomic regions under selection and the environmental parameters responsible for their selection (e.g., honeybees in the Iberian climate Henriques et al., 2018; dogs with the elevation Gou et al., 2014; or goats in the Moroccan climate Benjelloun et al., 2015). In the meantime, the context of global changes lead to an increasing interest on the ability of populations to develop short-term responses to face e.g., climate variations, as well as the interaction between short- and long-term mechanisms. Indeed, populations might quickly react to an environmental change through migration towards more favorable environments, which is well documented in the global warming context (Chen et al., 2011; McDonald et al., 2012). In the short-term, phenotypic plasticity can also increase the ability of organisms to cope with environmental changes. The ability of a genotype to produce, under different conditions, different phenotypes that are best fitted to the environment involves non-genetic mechanisms by which favorable variations can be acquired and even sometimes transmitted. They consist in epigenetic variations, parental effects, ecological and cultural variations (Danchin et al., 2011). One of such mechanism is epigenetics, notably molecular modifications of chromatin without modification of the DNA sequence itself (including e.g., DNA methylation or histone modification) which affects gene expression (Gibney and Nolan, 2010). Until now, most of the studies pointing out the role of molecular epigenetic mechanisms were performed in stress-controlled conditions on plant or animal models. In this case, the effect of one varying parameter on the physiological or epigenetic responses of the organisms is assessed. For example, the glyphosate herbicide injury on Arabidopsis thaliana (Kim et al., 2017), or the maternal diet in mice (Cooney et al., 2002) were shown to affect the methylome.

Since the end of the 2000s, the role of these epigenetic processes in ecology and evolution is increasingly being studied (Bossdorf et al., 2008). In this context, epigenetic marks were analyzed in natural populations. For example, the search for epigenetic differentiation by contrasting natural populations of Lilium bosniacum living in different habitat conditions (Zoldoš et al., 2018), or in wild baboons (Papio cynocephalus) populations with different food resources (natural fodder in a savanna environment or human food scraps, Lea et al., 2016).

Like wild species, farm animals will have to face future climate change, and understanding the genetic and epigenetic mechanisms responding to environmental changes, especially for free-ranging populations, is a key issue for the conservation of Farm Animal Genomic Resources (Bruford et al., 2015). Our study focuses on small domestic ruminants, namely sheep (Ovis aries) and goats (Capra hircus), which are a very good model for studying local adaptation. Indeed they spread all over the world for about 10,000 years from the domestication center in the Middle East, through human migration and commercial trades (Zeder, 2008). Until recent centuries, they were traditionally managed in a sustainable way leading to populations well adapted to the large variety of environments worldwide (Taberlet et al., 2008). Nowadays, the effect of the environment on farm animals strongly depends on breeding conditions. Intensively farmed animals are rather preserved from environmental variation because they are generally kept indoors with access to quality supply, in contrast to extensive livestock which are raised outdoors with little additional water and feed. Since the rise of the breed concept about 200 years ago, the selection for morphological and productivity traits improved the productivity, but resulted in the loss of genetic variability in many breeds, compromising their ability to adapt to future environmental changes (Taberlet et al., 2008). Thus, in the context of climate change gathering “epigenomic information to be integrated with phenotypic and genomic data to scrutinize the biological basis for adaptation and plasticity/resilience in livestock populations” is a key question in farm animal genomics research (Bruford et al., 2015). Until now, epigenetic analyses on livestock have been performed almost exclusively in stress-controlled conditions in relation with agronomic traits such as reproduction (Lan et al., 2013), production (Peñagaricano et al., 2014; X. Wang et al., 2015) and disease sensibility (Doherty et al., 2014). To our knowledge, only the study by Sevane et al. (2018) searched for DNA methylation in relation with climate variation by comparing cattle adapted to tropical and temperate climates.

In this context we aimed at characterizing both epigenetic and genetic variations with regards to climate variation in a single framework. Our study focused on sheep and goats sampled in Morocco where farming is mainly carried out in traditional ways with a strong impact of environmental conditions upon animals. Within each species, our goal was to identify candidate genomic regions involved in both genetic (selection on Single Nucleotide Polymorphism, SNP) and epigenetic (differential DNA methylation) mechanisms related to environmental variations, using whole genome comparisons of groups of individuals from locations with contrasted annual temperature range.

 Materials and Methods Sampling

Morocco has very contrasted climates including hot/cold desert and semi-arid, warm-summer Mediterranean, dry-summer subarctic (Born et al., 2008). In this context, we sampled sheep (Ovis aries) and goats (Capra hircus) from locations with contrasted annual temperature variations (Figure 1, bioclimatic variable BIO7 from the WorldClim dataset, www.worldclim.org/bioclim, Hijmans et al., 2005). This variable is the difference between the maximum temperature of the warmest month and the minimum temperature of the coldest month. It is correlated with several temperature variables, sunshine and elevation (Supplementary Table S1). Twenty-one unrelated goats (10 and 11 from low and high annual temperature ranges, respectively) and 22 unrelated sheep (12 and 10 from low and high annual temperature ranges, respectively) were selected (Figure 1; Supplementary Table S2) from the samples collected between January 2008 and March 2012 in the frame of the NextGen European project (Grant Agreement no. 244356), in accordance with ethical regulations of the European Union Directive 86/609/EEC, and for which Whole Genome Sequences were already available (projects.ensembl.org/nextgen/). The average annual temperature range between the two groups is 20 vs. 39.5°C, which is almost twice for the animals living in Eastern Morocco.

Variation of the annual temperature range and sampling of small ruminants in Morocco. Distribution of sheep and goats sampled at each end of the temperature gradient in Morocco. The color gradient corresponds to the variation of the annual temperature range (BIO7).

 Whole Genome Sequences Dataset

We extracted the states of SNPs from the variant files generated by the NextGen project. Whole Genome Sequences (WGS) at 12–14 X coverage corresponding to the studied individuals were retrieved from the European Nucleotide Archive (www.ebi.ac.uk/ena, accession code PRJEB7436). The sequences, produced and filtered following the protocol in Alberto et al. (2018), were aligned on the reference genomes CHIR v1.0 (GenBank assembly GCA_000317765.147, Dong et al., 2013) for goats and OAR v3.1 (GenBank assembly GCA_000298735.1, Jiang et al., 2014), for sheep.

 Detection of Selection Signatures

We kept SNPs with a missing data rate inferior to 30% and a minor allele frequency above 1%. Then we calculated the FST value (Weir and Cockerham, 1984) between the two groups (low vs. high temperature annual range for each species) for each SNP using the weir-fst-pop function in vcftools version 0.1.16 (Danecek et al., 2011). Candidate genomic regions were defined by 1) merging SNPs with the highest FST values (i.e., top 0.025% and top 0.001% for the pathway and gene approaches, respectively) closer than 5,000 bases from each other and then 2) adding 2,500 bases at each side of the region obtained. At the end 3) the genes overlapping these regions were extracted.

To annotate the genomic regions, two gene ontology analyses were performed with the GO Ontology database released on 2020-08-10 (Carbon and Mungall, 2018). The Bos taurus reference was used to annotate the candidate genes. We do not expect any major issue due to the use of B. taurus as a reference as the genome synteny between these three Bovidae species is very high (see Jiang et al., 2014).

1) Enrichment analysis: it was performed on the list of genes extracted with the 0.025% threshold. For that, we used the statistical over-representation test in PANTHER (version 15.0, Mi et al., 2019) and the Bos taurus reference. Each biological process with a significant Fisher’s exact test and a False Discovery Rate (FDR) lower than 0.05 was annotated. We kept the more basal significant biological pathway and searched the GO terms with the highest hierarchical level (i.e., child term of “biological process”) associated.

2) Candidate gene identification: we extracted the GO term corresponding to the top genes using a more stringent threshold (0.001%) to avoid the false positive and highlight the most impacted genes. In addition, a specific gene bibliography was done using Google Scholar with the gene’s name and “goat/sheep” or “livestock” or “mammals” as keywords, to determine what could be their effect on livestock phenotypes. These were grouped into large phenotype categories such as “Environmental perception,” “Immunity,” “Production” (genes involved in milk or meat production), “Reproduction” and “Other.”

DNA Methylated Regions DNA Extraction, Library Preparation and Sequencing

DNA was extracted from ear biopsies collected by the Nextgen consortium (see above) using the isolation of genomic DNA from tissues protocol from QIAamp® DNA micro kit (QIAgen, Germantown, MD, United States). Methylated DNA immunoprecipitation sequencing (MeDIP-Seq) was performed at the GeT-PlaGe core facility, INRAE Toulouse (http://www.get.genotoul.fr). Methylated DNA libraries were prepared according to Bioo Scientific’s protocol using the Bioo Scientific NEXTflex™ Methyl-seq Library Prep Kit for Illumina Sequencing. Briefly, DNA was fragmented by sonication on a covaris M220, size selection was performed using AMPure XP beads and adapters were ligated to be sequenced. 3 µg of each library was diluted in 450 µl of TE buffer, denatured in boiling water for 10 min and immediately cooled in ice for 10 min. Then 50 µL of 10X concentrated IP buffer was added to the mixture as well as 1 µg of anti-5-methylcytosine monoclonal antibody (clone 33D3, Diagenode®). After overnight incubation at 4°C with agitation, the DNA-antibody complexes were purified using 40 µL of “Dynabeads-ProteinG” previously washed in PBS. The DNA-antibody complexes were washed twice in 700ul of IP buffer and then the DNA was purified using the iPure kit (Diagenode®). Control or input DNA samples were not coupled to the antibodies and were therefore not enriched in epigenetic markers. Then, 12 cycles of PCR were performed. Library quality was assessed using an Advanced Analytical Fragment Analyzer and libraries were quantified by QPCR using the Kapa Library Quantification Kit. MeDIP-Seq was performed on an Illumina® HiSeq 2500 with the Illumina Reagent Kits. At the end, the sequenced reads correspond to fragment of methylated DNA.

 Dataset Preparation

Illumina paired-end reads (mean size of 125 bp) from all animals were aligned on the same reference genomes as WGS data by using the default parameters of BWA-MEM version 0.7.12 (Li and Durbin, 2009). The function multicov from bedtools version 2.29.2 (Quinlan and Hall, 2010) was used to report the MeDIP-Seq coverage of alignment files inside windows of 100 non-overlapping base pairs along the genome. The same function was used to extract the WGS coverage. In these windows we also extracted the number of CpG from the reference genomes by using the MEDIPS package (Lienhard et al., 2014) on R software version 3.5.3 (R Core Team, 2019).

We compared genome wide coverage profiles of MeDIP-Seq by using the MEDIPS.correlation function from the MEDIPS package with the uniq, extent, and shift parameters equal to 0 and the window size parameter equal to 100, to keep all MeDIP-Seq reads in each 100 bases window on the genome. We verified with a pairwise Pearson test that the number of read by windows for one animal had a correlation coefficient greater than 0.6 with all other individuals of the same species.

The numbers of WGS and MeDIP-Seq reads in each window were normalized between individuals with default parameters of calcNormFactors function (edgeR package, Robinson et al., 2010). Then, normalized MeDIP-Seq read counts were divided by normalized WGS read counts in order to get rid of the effect of sequencing coverage variations (e.g., due to repeated regions). Normalized values have been rounded to two decimals. Four outlying windows (2 for each species) with 10 times more MeDIP-Seq reads than the others were removed as this may be due to PCR or alignment errors.

 Whole Epigenome Correlations

For each species, correlations between the number of CpG and that of MeDIP-Seq reads found in the 100 bp window were tested using a Pearson’s correlation test. We used the same test to test the correlation between species for each of these two variables. For that we used the correspondence between Capra and Ovis orthologous regions resulting from the cross-alignment between the two reference genomes CHIR_1.0 and OAR v3.1 performed by Alberto et al. (2018). Results were represented as density plots performed with the hexbin package (Carr et al., 2019).

 Detection of Differentially Methylated Regions

We used the edgeR package to detect DMRs. Starting from non-standardized MeDIP-Seq read counts, we kept the 100 bp windows for which at least 10% of the samples had at least 1 count-per-million (1 read per million windows, this accommodates differences in library sizes between sample) and we normalized the read counts with the calcNormFactors function with default parameters. We estimated common dispersion with the estimateGLMRobustDisp function which mitigates the influence of outliers (X. Zhou et al., 2014) and we tested the occurrence of differential methylation between groups with the default parameters of the exactTest function. The p-values obtained were converted into q-values with the Bioconductor package qvalue (Storey et al., 2015) and the ones showing q-values inferior to 0.01 were selected for the rest of the analysis. Finally, we merged significant windows distant than less than 5,000 bases into regions.

In order to delineate the peaks of methylation occurring for some of the individuals within these regions, we used the changepoint package (Killick and Eckley, 2014) which detects breakpoints in time series. We looked for all possible breakpoints within a region using the cpt.meanvar function that identify change in mean and variance of the sum of coverages for all animals of each group. We kept the region between two breakpoints that contained at least one window detected as significant with edgeR, thus, the delimited regions were merged with the bedr package (Waggott et al., 2017) when they overlapped. Then, we used a Wilcoxon test to detect peaks with different methylation levels between the two groups. Among those peaks, we discarded those where the number of MeDIP-Seq reads was related to the DNA coverage in the WGS data, because here high MeDIP-Seq counts would result from the high DNA coverage (e.g., due to repeated sequences). The regions containing the remaining peaks were mapped on the reference genomes and included or overlapping genes were identified. We used the Functional classification viewed in gene list analysis in PANTHER to extract the Biological Pathways annotated in the Gene Ontology Database associated to these genes.

 Results Detection of Selection Signatures

The FST values at each genome position and their overall distribution are presented for goats and sheep in Supplementary Figures S1, S2, respectively. Almost 94,5% of the FST values are less or equal to 0.1 for both species. The top 0.025 and 0.001% thresholds correspond to FST value of 0.4 and 0.55 respectively.

Enrichment Analysis

The 5,000 SNPs for goat and 6,000 SNPs for sheep which harbored the highest 0.025% FST values (Supplementary Figures S3, S4) corresponded to 442 and 489 genes, respectively, putatively under selection.

The enriched pathways were associated with 8 and 10 GO terms in goat (Figure 2, light green color; Supplementary Table S3) and sheep (Figure 2, dark green color; Supplementary Table S4), respectively. Eight were common to both species.

GO terms putatively involved in the adaptation and acclimation of small ruminants to temperature variations. Number of genes or child GO terms per basal GO terms for candidate gene identification (under selection -oranges- and differentially methylated -yellow/brown-) and enrichment analysis (greens), respectively.

Among the 0.025% top genes (i.e., 931 genes selected for the enrichment analysis), 41 were orthologous between sheep and goats. When we did the same enrichment analysis (described above) with the 41 genes shared by sheep and goats we obtained 21 GO terms. Four of them had the “Nervous system development” in their parental terms; and the last three were associated with cellular processes, one of which was related to synapses organization.

 Candidate Gene Identification

The top 0.001% FST values corresponded to 194 SNPs for goat and 146 SNPs for sheep distributed in 56 and 42 genomic regions, respectively. Their location in the species genome is detailed in Supplementary Tables S5, S6 for goat and sheep, respectively. In goats, 34 of the 56 regions were intergenic and the others were related to 22 genes (Table 1). For sheep, we obtained 24 intergenic regions and 18 genes (Table 1), which differed between the two environmental groups. No genes were common between the two species.

Phenotype categories in which the candidate genes (top 0.001%) are involved.

Species Phenotype categories Genes
Goat Environmental perception TMTC2, EDIL3, SASH1
Immunity —
Production CP, SLC9A9, PPFIA2, TMTC2, GALNTL6, CTSB, SASH1, PLCG2, MYADM
Reproduction ADGRB3, USH1C, PLCG2, MYADM, SERPINB7
Other TRABD2B, MVB12B, KLF12, TNFSF9, CALHM3, SORCS1, RAB30
Unknown LOC102184299, LOC106503718
Sheep Environmental perception NOX3, KSR2
Immunity SAMD12, SEMA5A, SEZ6L
Production BMPER, KSR2, PTPRE
Reproduction GCSAML, KSR2, BMPER
Other AGO3, IQCJ, KCNG3, ANO6, CASP8AP2, CNKSR3, ENGASE, KCNA4
Unknown LOC101107868, LOC101119001

The 22 genes were linked to 15 GO terms in goat (Figure 2, light orange color; Supplementary Table S5) and the 18 genes in sheep were associated with 14 GO terms (Figure 2, dark orange color; Supplementary Table S6). Twelve of these GO terms were common to both species.

In livestock, several of these genes are known to be involved in phenotypic changes in relation with the perception of the environment (3 genes in goat, 2 in sheep), immunity (1 in goat, 3 in sheep), production (9 in goat, 2 in sheep), reproduction (5 in goat and 3 in sheep). The “Other” category brings together the genes whose effect on the phenotype was not documented in livestock species. We also detected genes that were uncharacterized (Table 1, Supplementary Tables S5, S6).

 Distribution of CpG and Methylated Regions in Both Species

The majority of the 22 million windows (i.e., 96% for each species) carried one methylated read, 3 and 2% more than one, and 1 and 2% none, for goat and sheep, respectively. Concerning the number of CpG, 47% of the windows did not have any, 50% had between 1 and 4, and 3% had more than 4 in both species.

Genomic windows of 100 bp that were orthologous between Capra and Ovis showed a correlation between species for the number of methylated read (R = 0.81, p-value < 2.10-16, Figure 3A) and the number of CpGs (R = 0.67, p-value < 2.10-16, Figure 3B) A large number of windows had few CpG and methylated reads (25% carried less than 3 CpGs and 1 methylated read). Within species, we found a correlation between the number of CpGs and that of methylated reads, R = 0.49 (p-value < 2.10-16, Figure 3C) for goats and R = 0.44 (p-value < 2.10-16, Figure 3D) for sheep.

Characterization and comparison of sheep and goats methylomes. Relationship between the number of (A). CpG, or (B). MeDIP-Seq reads per 100 bp window between the two species across the genome. Relation between the number of MeDIP-Seq reads and the number of CpG per 100 bp widow for (C). goat, and (D). sheep. The colors represent the number of points and the black line (y = x) is represented.

For the detection of DMRs, 20,256 windows in goat and 16,281 windows in sheep were kept out of more than 24 millions. The edgeR analysis revealed 20 and 7 windows with significant differences in methylation levels between the two groups (i.e., low vs high annual temperature variations) in goat and sheep, respectively. Merging candidate windows closer than 5,000 bp, we obtained 8 and 4 regions in which we detected 8 and 5 methylation peaks for goat and sheep, respectively. From these, 4 and 2 peaks respectively showed significant differential MeDIP-Seq coverage between the 2 groups (low vs. high temperature annual range). For 4 peaks, the differential MeDIP-Seq coverage was concomitant to a differential WGS coverage, indicating that the differentiation would probably result from repetition polymorphism. Then only one region per species was kept as resulting from a differential methylation. They were associated to AGPTA4 in goat and SLIT3 in sheep (Figure 4).

Differentially methylated regions putatively involved in temperature acclimation in small ruminants. DMRs between genomes from locations with low (blue) and high (red) temperature annual range for (A). goat or (B). sheep. Each line represents the variations of the level of DNA methylation along a gene for one individual; the grey bar represents part of the gene; the two black lines delineate the peak on which the differential level of methylation was tested. Boxplots represent the level of methylation (MeDIP-Seq coverage) for each group compared with a Wilcoxon test (Significance codes: ** = p-value ≤ 0.01, * = p-value ≤ 0.05).

The GO terms associated with these genes were “Cellular process,” and “Metabolic process” for AGPAT4 (Figure 2, yellow color); and “Biological regulation,” “Cellular process,” “Developmental process,” “Growth,” “Locomotion,” “Metabolic process,” “Multicellular organismal process,” “Reproduction,” “Reproductive process,” “Response to stimulus,” “Rhythmic process,” and “Signaling” for SLIT3 (Figure 2, brown color).

 Discussion

To our knowledge, this is the first study to investigate simultaneously the genetic and epigenetic variations putatively involved in the adjustment of phenotypes to field conditions in small ruminants. Until now, few publications jointly addressed both mechanisms. Some studies highlighted significant association between local genetic variation and the presence of DMRs like that of Eichten et al. (2013) who found that 51% of the DMRs between maize populations were associated with local SNPs. Foust et al. (2016) and Herrera and Bazaga (2010) found a correlation between AFLP and MS-AFLP marks in two salt marsh perennials (Spartina alterniflora and Borrichia frutescens) and violet populations (Viola cazorlensis), respectively. In our study, we investigated both genetic and epigenetic variations in sheep and goats originating from field populations. The sampling was carried out in order to avoid a breed effect as much as possible. However, given the specific distribution and adaptation of some breeds, several Ghazalia goats and D’man sheep were part of the groups submitted to high temperature annual range (c.a. 39.5°C) as opposed to that with low temperature annual range (c.a. 20°C). They represented less than half of the individuals (3/11 and 4/10 in goats and sheep, respectively) and had similar methylation and FST patterns to the other breeds from the same group. Thus, the mechanisms inferred would be shared by the variety of Moroccan breeds represented in our sampling. Our study points out several genes which polymorphism or methylation level is contrasted between the two groups and may contribute to the adaptation and acclimatization to these specific environments. The candidate genes putatively under selection are primarily involved in environmental perception, immunity, production and reproduction, and their variations could be driven by different environmental factors. First, a direct effect of high temperatures might operate through thermoregulation by reducing body temperature, through the decrease of metabolic heat and the increase of heat dissipation. This might be obtained by avoiding solar radiations (Al-Tamimi, 2007), increasing respiration or sweating and reducing food intake (Marai et al., 2007; Sejian et al., 2010). These adjustments might reduce the energy allocated to other biological processes including productivity (Tao et al., 2011; Dado-Senn et al., 2019) or reproduction (Meyerhoeffer et al., 1985; Monteiro et al., 2016; Krishnan et al., 2017). In relation with temperature annual range, we found candidate genes involved in production, including milk production for CP, MYAMD, PTPRE, SASH1 and SLC9A9 (Hussein and Staufenbiel, 2012; Dong et al., 2015; Laodim et al., 2017; Shin et al., 2017; Xiang et al., 2017), mammary gland development for PPFIA2 (Mota et al., 2020) and morphology for TMTC2 (Abo-Ismail et al., 2017; Z.-H. Fang and Pausch, 2019). In addition, CTSB, GALNTL6, and KSR2 genes were related to meat production (Russo et al., 2002; Doran et al., 2014; Puig-Oliveras et al., 2014; Lukić et al., 2020), and BMPER and KSR2 with body size (Puig-Oliveras et al., 2014; Zhao et al., 2015; Lukić et al., 2020). Genes involved in reproduction were also impacted. ADGRB3, GCSAML, and KSR2 affect fertility (Pimentel et al., 2011; X. Wang et al., 2015; Mohammadi et al., 2020), while PLCG2 and SERPINB7 affect gestation (Vallée et al., 2003; Ponsuksili et al., 2012). Moreover offspring survival can also be impacted through litter size via BMPER and USH1C (Kwon et al., 2016; Bolormaa et al., 2017; Mota et al., 2020) and offspring weight with MYADM (Gonzalez et al., 2020). Also, we found two candidate genes related to environmental perception. NOX3 is involved in the response to stresses in general (Lambeth, 2004) and has been related with heat stress in chicken (Aggrey et al., 2018), and KSR2 plays a role in the energy balance by adapting feeding behavior and thermogenesis in mice (Guo et al., 2017). Second, environmental parameters covarying with the annual temperature range (BIO7) may be prominent in driving the selection on the genes which polymorphism was contrasted between the two groups. TMTC2 and EDIL3 are involved in adaptation to high altitude in sheep and cattle (Ji Yang et al., 2016; Pierce et al., 2020), this later variable being correlated with BIO7 in Morocco. In addition, BIO7 variable is correlated with sunshine which impacts directly skin pigmentation, which would explain the detection of SASH1, which regulates melanocytes migration in relation with the production and transfer of melanin to protect skin from UV radiations (D. Zhou et al., 2013). Third, the temperature regime may have indirect impacts on phenotypes by driving other environmental factors. Especially, it determines plant resources (i.e., food quality and quantity) and is related to water availability, and thus strongly impacts the diet of free ranging livestock. The lack of resources has a negative impact on, e.g., reproduction (Martin et al., 2004; J. J. Robinson et al., 2006) and lactation (Razzaghi et al., 2016). In contrast, supplementation, which has a positive effect on milk yield and composition in sheep, could temperate this negative impact on lactation (Caroprese et al., 2011). As a candidate gene, we found PLCG2 which is regulated in the cattle liver according to feed intake (Salleh et al., 2017). Likewise, genes involved in immune functions may be selected by multi-factorial drivers related to nutrition or the pathogenic context in relation with climate (Bett et al., 2017). This could be the case for SAMD12 playing a role in pancreatitis in human (Giri et al., 2016), SEZ6L associated with foot-and-mouth disease (Lee et al., 2015) or SEMA5A with mastitis in cattle (Sugimoto et al., 2006).

Besides the candidate genes putatively involved in local adaptation through allelic selection, we identified two strong candidates in relation with an epigenetic response. The global epigenomic context was similar in sheep and goats, and the level of DNA methylation along the genome was correlated to the CpG content. This is consistent with what is known about mammalian epigenomes where the CpG methylation generally covers the entire nuclear genome, with the exception of CpG rich regions near the promoters of active genes (Suzuki and Bird, 2008; Feng et al., 2010). We found two genes differentially methylated according to the provenance of the individual, namely AGPAT4 (1-acylglycerol-3-phosphate O-acyltransferase 4) in goats, and SLIT3 (Slit Guidance Ligand 3) in sheep. AGPAT4 is part of the 1-acyl-sn-glycerol 3-phosphate acyltransferases family which is involved in triglyceride synthesis (Lu et al., 2005; Takeuchi and Reue, 2009), with a role in the milk composition of ruminants (Bionaz and Loor, 2008; He et al., 2011; Jie Yang et al., 2016). It has an organ-related level of expression in human (Lu et al., 2005; Prasad et al., 2011) and plays a role in muscle development and meat quality in cattle (X. Fang et al., 2017). In addition, the variation of AGPAT genes expression in the mammary gland has been related to milk fat depression in dairy cows (Peterson et al., 2003). On its side, SLIT3 takes place into the SLIT-ROBO interactions which are involved in the development of the nervous system (K. H. Wang et al., 1999; Andrews et al., 2007), cardiac development (Liu et al., 2003; Mommersteeg et al., 2015) and reproduction, especially menstrual cycle (Dickinson and Duncan, 2010) and fertility (Amundson et al., 2015). SLIT3 was also shown to be differentially methylated and upregulated in goat mammary glands during lactation compared to the dry period (Zhang et al., 2017). Such DMRs related to milk quality and lactation might result directly from heat stress. Indeed, it was shown that heifers submitted to a thermic shock during fetal development had smaller mammary gland alveoli, resulting in lower milk production (Skibiel et al., 2018; Dahl et al., 2019). The DMRs observed could also result from indirect effects of the climate through a variation of the diet resulting from a shift in plant resources. Peterson et al. (2003) showed that the concentrate forage diet influences the mRNA abundance of AGPAT genes in mammary glands in dairy cows, involving a modification of milk fat. The regulation of the SLIT-ROBO pathway is also affected by the diet (Amundson et al., 2015). Indeed, the diet can directly influence the levels of DNA methylation (Anderson et al., 2012) either due to direct exposure of the individuals or to parental effects. For example, a heat-stress-induced perturbation of mammary glands development in the fetus, related to differential methylation and differential gene expression (Skibiel et al., 2018), results in effects persisting over the lactations and in the next generation too (Ouellet et al., 2020). Maternal nutrition also affects DNA methylation and gene expression of lambs, especially in muscle and adipose tissues (Lan et al., 2013; Peñagaricano et al., 2014; Namous et al., 2018). Obviously, we found differential methylation in genes related to lactation in ear biopsies, not mammary glands. However, sweat glands localized in the dermis and mammary glands have a common tissue origin, and recent studies showed that the most similar tissues had comparable CpG methylation patterns (Lokk et al., 2014). Given the epigenetic memory, embryonic patterns may be retained in adult tissues (Hon et al., 2013). In this case, an early stress (e.g., in utero) could modify the methylation pattern of embryonic cells, which could be retained in structures with a common origin. Further characterization of the methylation patterns in different organs, including mammary glands, is needed to confirm the possible role of these mechanisms.

We found that candidate genes for local adaptation in response to different temperature regimes were related to a large variety of biological pathways which were mainly involved in functions such as immunity, perception of the environment and reproduction, or related to productivity traits. While these functions were affected in both sheep and goats, the best candidate genes (i.e., top 0.001%) identified were different, confirming that alternative genes might be selected under the same environmental constraints for close species (e.g., Benjelloun, 2015; Raeymaekers et al., 2017) or even populations of the same species (Manceau et al., 2010; Benjelloun et al., 2015). However, out of the 0.025% top genes (i.e., 931 genes selected for the enrichment analyses), 41 were orthologous between sheep and goats and mainly related to the nervous system development. Interestingly, the present study also investigated epigenetic variations in the individuals displaying candidate genes for local adaptation. Our deliberately stringent approach to avoid false positive, which identified only one candidate DMR per species, prevents a global view but already points out similarities between traits impacted by genetic and epigenetic mechanisms. In sheep, genes impacted by selection such as BMPER, GCSAML, KSR2, PTPRE affect embryo quality, litter size and lactation, which are also influenced by SLIT3 that was differentially methylated. Likewise, in goats, AGPAT4 that was differentially methylated has an effect on lipid metabolism, meat quality and the presence of lipids in milk, these traits being influenced by several candidate genes such as CTSB, MYADM, SERPINB7, USH1C, GALNTL6 and CALHM3.

 Conclusion

This study showed differentially methylated genes and candidate genes under selection related to a differential regime of temperatures underwent by sheep and goats. Genetic and epigenetic mechanisms share several target functions including lactation and milk composition. These adjustments of phenotypes to the two contrasted environments may result from drivers related to temperature either directly or indirectly through e.g., the variation of plant resources conditioning the diet. We do not yet have enough evidence to say that these results reflect the joint effects of selection and epigenetic regulation on the same traits. However, such combined effects might play an important role in evolution (Richards, 2011), and we can wonder whether DMRs might cushion phenotypic variations due to prior selection of genes or reflect the persistence of adaptive plasticity as long as adaptation is not complete? Deciphering the relative impact of climate and e.g., differential diets on the genetic and epigenetically-induced variations of phenotypes will be a major stake for developing future breeding strategies in the context of global change.

 Data Availability Statement

The genetic dataset (fastq of Whole Genome Sequencing) analyzed in this study can be found in the European Nucleotide Archive (www.ebi.ac.uk/ena, accession code PRJEB7436). The epigenetic dataset (fastq of the MeDIP-seq) generated for this study can be found in the European Nucleotide Archive (www.ebi.ac.uk/ena, accession code PRJEB47531).

 Ethics Statement

The animal study was reviewed and approved by European Union Directive 86/609/EEC.

 Author Contributions

FP and FB supervised the project. FP and BB contributed to the design project. BB conducted the sampling of Moroccan small ruminants. FA, FP, and LD performed DNA extraction and managed the DNA samples. LD, CG, MK, PD, and FB performed the analyzes. LD and FP wrote the paper.

 Funding

We thank the Ministère de l’Enseignement Supérieur, de la Recherche et de l’Innovation for financing LD. The LECA is part of the Labex OSUG@2020 (ANR10 LABX56). This study was partly funded by the projects Climgen (Program FACE ERA-NET PLUS 2014, ANR- 14-JAFC-0002) and MethylGoat (AGIR-2014, Univ. Grenoble Alpes). Computations were performed using the CIMENT infrastructure, which is supported by the Auvergne-Rhône-Alpes region (CPER07_13 CIRA).

 Conflict of Interest

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

 Publisher’s Note

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

We thank Florence Plisson and Sylvie Veyrenc for their help during the DNA manipulations (extraction, MeDIP-Seq); and Hervé Acloque for help in setting up the MeDIP-Seq approach. This work was performed in collaboration with the GeT core facility, Toulouse, France (http://get.genotoul.fr), and was supported by France Génomique National infrastructure, funded as part of “Investissement d’avenir” program managed by Agence Nationale pour la Recherche (contract ANR-10-INBS-09). All the computations presented in this paper were performed using the GRICAD infrastructure (https://gricad.univ-grenoble-alpes.fr), which is partly supported by the Equip@Meso project (reference ANR-10-EQPX-29-01) of the programme Investissements d’Avenir supervised by the Agence Nationale pour la Recherche.

 Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2021.745284/full#supplementary-material

References Abo-Ismail M. K. Brito L. F. Miller S. P. Sargolzaei M. Grossi D. A. Moore S. S. (2017). Genome-wide Association Studies and Genomic Prediction of Breeding Values for Calving Performance and Body Conformation Traits in Holstein Cattle. Genet. Sel Evol. 49 (1), 82. 10.1186/s12711-017-0356-8 Aggrey S. E. Habashy W. S. Milfort M. C. Adomako K. Fuller A. L. Rekaya R. (2018). “Molecular and Cellular Mechanisms that Underlie Genes and Antioxidant Enzyme Activities in Meat-type Birds during Heat Stress,” in World Congress on Genetics Applied to Livestock Production, Auckland, New Zealand. Al-Tamimi H. J. (2007). Thermoregulatory Response of Goat Kids Subjected to Heat Stress. Small Ruminant Res. 71 (1), 280–285. 10.1016/j.smallrumres.2006.04.013 Alberto F. J. Boyer F. Orozco-terWengel P. Streeter I. Servin B. de Villemereuil P. (2018). Convergent Genomic Signatures of Domestication in Sheep and Goats. Nat. Commun. 9 (1), 813. 10.1038/s41467-018-03206-y Amundson O. L. Fountain T. H. Larimore E. L. Richardson B. N. McNeel A. K. Wright E. C. (2015). Postweaning Nutritional Programming of Ovarian Development in Beef Heifers1,2. J. Anim. Sci. 93 (11), 5232–5239. 10.2527/jas.2015-9067 Anderson O. S. Sant K. E. Dolinoy D. C. (2012). Nutrition and Epigenetics: An Interplay of Dietary Methyl Donors, One-Carbon Metabolism and DNA Methylation. J. Nutr. Biochem. 23 (8), 853–859. 10.1016/j.jnutbio.2012.03.003 Andrews W. D. Barber M. Parnavelas J. G. (2007). Slit?Robo Interactions during Cortical Development. J. Anat. 211 (2), 188–198. 10.1111/j.1469-7580.2007.00750.x Benjelloun B. Alberto F. J. Streeter I. Boyer F. d. r. Coissac E. Stucki S. (2015). Characterizing Neutral Genomic Diversity and Selection Signatures in Indigenous Populations of Moroccan Goats (Capra hircus) Using WGS Data. Front. Genet. 6, 107. 10.3389/fgene.2015.00107 Benjelloun B. (2015). “Diversité des génomes et adaptation locale des petits ruminants d’un pays méditerranéen: Le Maroc,”. doctoral thesis (Grenoble (France): Université Grenoble Alpes). http://www.theses.fr/2015GREAV011. Bett B. Kiunga P. Gachohi J. Sindato C. Mbotha D. Robinson T. (2017). Effects of Climate Change on the Occurrence and Distribution of Livestock Diseases. Prev. Vet. Med. 137, 119–129. 10.1016/j.prevetmed.2016.11.019 Bionaz M. Loor J. J. (2008). ACSL1, AGPAT6, FABP3, LPIN1, and SLC27A6 Are the Most Abundant Isoforms in Bovine Mammary Tissue and Their Expression Is Affected by Stage of Lactation. J. Nutr. 138 (6), 1019–1024. 10.1093/jn/138.6.1019 Bolormaa S. Brown D. J. Swan A. A. van der Werf J. H. J. Hayes B. J. Daetwyler H. D. (2017). Genomic Prediction of Reproduction Traits for Merino Sheep. Anim. Genet. 48 (3), 338–348. 10.1111/age.12541 Born K. Christoph M. Fink A. H. Knippertz P. Paeth H. Speth P. (2008). “Moroccan Climate in the Present and Future: Combined View from Observational Data and Regional Climate Scenarios,” in Climatic Changes and Water Resources in the Middle East and North Africa. Editors Zereini F. Hötzl H. (Springer Berlin Heidelberg), 29–45. 10.1007/978-3-540-85047-2_4 Bossdorf O. Richards C. L. Pigliucci M. (2008). Epigenetics for Ecologists. Ecol. Lett. 11, 106–115. 10.1111/j.1461-0248.2007.01130.x Bruford M. W. Ginja C. Hoffmann I. Joost S. Orozco-terWengel P. Alberto F. J. (2015). Prospects and Challenges for the Conservation of Farm Animal Genomic Resources, 2015-2025. Front. Genet. 6, 314. 10.3389/fgene.2015.00314 Carbon S. Mungall C. (2018). Gene Ontology Data Archive [Data Set]. Zenodo. 10.5281/zenodo.3980761 Caroprese M. Albenzio M. Bruno A. Fedele V. Santillo A. Sevi A. (2011). Effect of Solar Radiation and Flaxseed Supplementation on Milk Production and Fatty Acid Profile of Lactating Ewes under High Ambient Temperature. J. Dairy Sci. 94 (8), 3856–3867. 10.3168/jds.2010-4067 Carr D. Lewin-Koh N. Maechler M. Sarkar D. (2019). Hexbin: Hexagonal Binning Routines. (R package version 1.28.0) [Computer software]. https://CRAN.R-project.org/package=hexbin. Chen I.-C. Hill J. K. Ohlemüller R. Roy D. B. Thomas C. D. (2011). Rapid Range Shifts of Species Associated with High Levels of Climate Warming. Science 333 (6045), 1024–1026. 10.1126/science.1206432 Cooney C. A. Dave A. A. Wolff G. L. (2002). Maternal Methyl Supplements in Mice Affect Epigenetic Variation and DNA Methylation of Offspring. J. Nutr. 132 (8), 2393S–2400S. 10.1093/jn/132.8.2393S Dado-Senn B. Skibiel A. L. Fabris T. F. Dahl G. E. Laporta J. (2019). Dry Period Heat Stress Induces Microstructural Changes in the Lactating Mammary Gland. PLOS ONE 14 (9), e0222120. 10.1371/journal.pone.0222120 Dahl G. E. Skibiel A. L. Laporta J. (2019). In Utero Heat Stress Programs Reduced Performance and Health in Calves. Vet. Clin. North America: Food Anim. Pract. 35 (2), 343–353. 10.1016/j.cvfa.2019.02.005 Danchin É. Charmantier A. Champagne F. A. Mesoudi A. Pujol B. Blanchet S. (2011). Beyond DNA: Integrating Inclusive Inheritance into an Extended Theory of Evolution. Nat. Rev. Genet. 12 (7), 475–486. 10.1038/nrg3028 Danecek P. Auton A. Abecasis G. Albers C. A. Banks E. DePristo M. A. (2011). The Variant Call Format and VCFtools. Bioinformatics 27 (15), 2156–2158. 10.1093/bioinformatics/btr330 Dickinson R. E. Duncan W. C. (2010). The SLIT-ROBO Pathway: a Regulator of Cell Function with Implications for the Reproductive System. Reproduction 139 (4), 697–704. 10.1530/REP-10-0017 Doherty R. Farrelly C. O. Meade K. G. (2014). Comparative Epigenetics: Relevance to the Regulation of Production and Health Traits in Cattle. Anim. Genet. 45, 3–14. 10.1111/age.12140 Dong Y. Xie M. Jiang Y. Xiao N. Du X. Zhang W. (2013). Sequencing and Automated Whole-Genome Optical Mapping of the Genome of a Domestic Goat (Capra hircus). Nat. Biotechnol. 31 (2), 135–141. 10.1038/nbt.2478 Dong Y. Zhang X. Xie M. Arefnezhad B. Wang Z. Wang W. (2015). Reference Genome of Wild Goat (capra Aegagrus) and Sequencing of Goat Breeds Provide Insight into Genic Basis of Goat Domestication. BMC Genomics 16 (431). 10.1186/s12864-015-1606-1 Doran A. G. Berry D. P. Creevey C. J. (2014). Whole Genome Association Study Identifies Regions of the Bovine Genome and Biological Pathways Involved in Carcass Trait Performance in Holstein-Friesian Cattle. BMC Genomics 15 (1), 837. 10.1186/1471-2164-15-837 Eichten S. R. Briskine R. Song J. Li Q. Swanson-Wagner R. Hermanson P. J. (2013). Epigenetic and Genetic Influences on DNA Methylation Variation in Maize Populations. Plant Cell 25 (8), 2783–2797. 10.1105/tpc.113.114793 Fang X. Zhao Z. Yu H. Li G. Jiang P. Yang Y. (2017). Comparative Genome-wide Methylation Analysis of Longissimus Dorsi Muscles between Japanese Black (Wagyu) and Chinese Red Steppes Cattle. PLoS ONE 12 (8), e0182492. 10.1371/journal.pone.0182492 Fang Z.-H. Pausch H. (2019). Multi-trait Meta-Analyses Reveal 25 Quantitative Trait Loci for Economically Important Traits in Brown Swiss Cattle. BMC Genomics 20 (1), 695. 10.1186/s12864-019-6066-6 Feng S. Cokus S. J. Zhang X. Chen P.-Y. Bostick M. Goll M. G. (2010). Conservation and Divergence of Methylation Patterning in Plants and Animals. Proc. Natl. Acad. Sci. 107 (19), 8689–8694. 10.1073/pnas.1002720107 Foust C. M. Preite V. Schrey A. W. Alvarez M. Robertson M. H. Verhoeven K. J. F. (2016). Genetic and Epigenetic Differences Associated with Environmental Gradients in Replicate Populations of Two Salt Marsh Perennials. Mol. Ecol. 25 (8), 1639–1652. 10.1111/mec.13522 Gibney E. R. Nolan C. M. (2010). Epigenetics and Gene Expression. Heredity 105 (1), 4–13. 10.1038/hdy.2010.54 Giri A. K. Midha S. Banerjee P. Agrawal A. Mehdi S. J. Dhingra R. (2016). Common Variants in CLDN2 and MORC4 Genes Confer Disease Susceptibility in Patients with Chronic Pancreatitis. PLOS ONE 11 (1), e0147345. 10.1371/journal.pone.0147345 Gonzalez M. Villa R. Villa C. Gonzalez V. Montano M. Medina G. (2020). Inspection of Real and Imputed Genotypes Reveled 76 SNPs Associated to Rear Udder Height in Holstein Cattle. J. Adv. Vet. Anim. Res. 7 (2), 234. 10.5455/javar.2020.g415 Gou X. Wang Z. Li N. Qiu F. Xu Z. Yan D. (2014). Whole-genome Sequencing of Six Dog Breeds from Continuous Altitudes Reveals Adaptation to High-Altitude Hypoxia. Genome Res. 24, 1308–1315. 10.1101/gr.171876.113 Guo L. Costanzo-Garvey D. L. Smith D. R. Neilsen B. K. MacDonald R. G. Lewis R. E. (2017). Kinase Suppressor of Ras 2 (KSR2) Expression in the Brain Regulates Energy Balance and Glucose Homeostasis. Mol. Metab. 6 (2), 194–205. 10.1016/j.molmet.2016.12.004 He C. Wang C. Chang Z. H. Guo B. L. Li R. Yue X. P. (2011). AGPAT6 Polymorphism and its Association with Milk Traits of Dairy Goats. Genet. Mol. Res. 10 (4), 2747–2756. 10.4238/2011.November.4.8 Henriques D. Wallberg A. Chávez-Galarza J. Johnston J. S. Webster M. T. Pinto M. A. (2018). Whole Genome SNP-Associated Signatures of Local Adaptation in Honeybees of the Iberian Peninsula. Sci. Rep. 8, 11145. 10.1038/s41598-018-29469-5 Herrera C. M. Bazaga P. (2010). Epigenetic Differentiation and Relationship to Adaptive Genetic Divergence in Discrete Populations of the Violet Viola Cazorlensis. New Phytol. 187 (3), 867–876. 10.1111/j.1469-8137.2010.03298.x Hijmans R. J. Cameron S. E. Parra J. L. Jones P. G. Jarvis A. (2005). Very High Resolution Interpolated Climate Surfaces for Global Land Areas. Int. J. Climatol. 25 (15), 1965–1978. 10.1002/joc.1276 Hon G. C. Rajagopal N. Shen Y. McCleary D. F. Yue F. Dang M. D. (2013). Epigenetic Memory at Embryonic Enhancers Identified in DNA Methylation Maps from Adult Mouse Tissues. Nat. Genet. 45, 1198–1206. 10.1038/ng.2746 Hussein H. A. Staufenbiel R. (2012). Variations in Copper Concentration and Ceruloplasmin Activity of Dairy Cows in Relation to Lactation Stages with Regard to Ceruloplasmin to Copper Ratios. Biol. Trace Elem. Res. 146 (1), 47–52. 10.1007/s12011-011-9226-3 Jiang Y. Xie M. Chen W. Talbot R. Maddox J. F. Faraut T. (2014). The Sheep Genome Illuminates Biology of the Rumen and Lipid Metabolism. Science 344 (6188), 1168–1173. 10.1126/science.1252806 Killick R. Eckley I. A. (2014). Changepoint: AnRPackage for Changepoint Analysis. J. Stat. Soft. 58 (3), 1–19. 10.18637/jss.v058.i03 Kim G. Clarke C. R. Larose H. Tran H. T. Haak D. C. Zhang L. (2017). Herbicide Injury Induces DNA Methylome Alterations in Arabidopsis. PeerJ 5, e3560. 10.7717/peerj.3560 Krishnan G. Bagath M. Pragna P. Vidya M. K. Aleena J. Archana P. R. (2017). “Mitigation of the Heat Stress Impact in Livestock Reproduction,” In Payan Carreira, R. (IntechOpen). Theriogenology Available at: https://www.intechopen.com/chapters/55491 . Kwon S. G. Hwang J. H. Park D. H. Kim T. W. Kang D. G. Kang K. H. (2016). Identification of Differentially Expressed Genes Associated with Litter Size in Berkshire Pig Placenta. PLOS ONE 11 (4), e0153311. 10.1371/journal.pone.0153311 Lambeth J. D. (2004). NOX Enzymes and the Biology of Reactive Oxygen. Nat. Rev. Immunol. 4 (3), 181–189. 10.1038/nri1312 Lan X. Cretney E. C. Kropp J. Khateeb K. Berg M. A. Peñagaricano F. (2013). Maternal Diet during Pregnancy Induces Gene Expression and DNA Methylation Changes in Fetal Tissues in Sheep. Front. Genet. 4, 49. 10.3389/fgene.2013.00049 Laodim T. Elzo M. A. Koonawootrittriron S. Suwanasopee T. Jattawa D. (2017). Identification of SNP Markers Associated with Milk and Fat Yields in Multibreed Dairy Cattle Using Two Genetic Group Structures. Livestock Sci. 206, 95–104. 10.1016/j.livsci.2017.10.015 Lea A. J. Altmann J. Alberts S. C. Tung J. (2016). Resource Base Influences Genome‐wide DNA Methylation Levels in Wild Baboons (Papio cynocephalus). Mol. Ecol. 25 (8), 1681–1696. 10.1111/mec.13436 Lee B.-Y. Lee K.-N. Lee T. Park J.-H. Kim S.-M. Lee H.-S. (2015). Bovine Genome-wide Association Study for Genetic Elements to Resist the Infection of Foot-And-Mouth Disease in the Field. Asian Australas. J. Anim. Sci. 28 (2), 166–170. 10.5713/ajas.14.0383 Li H. Durbin R. (2009). Fast and Accurate Short Read Alignment with Burrows-Wheeler Transform. Bioinformatics 25 (14), 1754–1760. 10.1093/bioinformatics/btp324 Lienhard M. Grimm C. Morkel M. Herwig R. Chavez L. (2014). MEDIPS: Genome-wide Differential Coverage Analysis of Sequencing Data Derived from DNA Enrichment Experiments. Bioinformatics 30 (2), 284–286. 10.1093/bioinformatics/btt650 Liu J. Zhang L. Wang D. Shen H. Jiang M. Mei P. (2003). Congenital Diaphragmatic Hernia, Kidney Agenesis and Cardiac Defects Associated with Slit3-Deficiency in Mice. Mech. Dev. 120 (9), 1059–1070. 10.1016/S0925-4773(03)00161-8 Lokk K. Modhukur V. Rajashekar B. Märtens K. Mägi R. Kolde R. (2014). DNA Methylome Profiling of Human Tissues Identifies Global and Tissue-specific Methylation Patterns. Genome Biol. 15 (4), 3248. 10.1186/gb-2014-15-4-r54 Lu B. Jiang Y. J. Man M. Q. Brown B. Elias P. M. Feingold K. R. (2005). Expression and Regulation of 1-Acyl-Sn-Glycerol- 3-phosphate Acyltransferases in the Epidermis. J. Lipid Res. 46 (11), 2448–2457. 10.1194/jlr.M500258-JLR200 Lukić B. Ferenčaković M. Šalamon D. Čačić M. Orehovački V. Iacolina L. (2020). Conservation Genomic Analysis of the Croatian Indigenous Black Slavonian and Turopolje Pig Breeds. Front. Genet. 11, 261. 10.3389/fgene.2020.00261 Manceau M. Domingues V. S. Linnen C. R. Rosenblum E. B. Hoekstra H. E. (2010). Convergence in Pigmentation at Multiple Levels: Mutations, Genes and Function. Phil. Trans. R. Soc. B 365 (1552), 2439–2450. 10.1098/rstb.2010.0104 Marai I. F. M. El-Darawany A. A. Fadiel A. Abdel-Hafez M. A. M. (2007). Physiological Traits as Affected by Heat Stress in Sheep-A Review. Small Ruminant Res. 71 (1–3), 1–12. 10.1016/j.smallrumres.2006.10.003 Martin G. B. Rodger J. Blache D. (2004). Nutritional and Environmental Effects on Reproduction in Small Ruminants. Reprod. Fertil. Dev. 16 (4), 491. 10.1071/RD04035 McDonald K. W. McClure C. J. W. Rolek B. W. Hill G. E. (2012). Diversity of Birds in Eastern North America Shifts north with Global Warming. Ecol. Evol. 2 (12), 3052–3060. 10.1002/ece3.410 Meyerhoeffer D. C. Wettemann R. P. Coleman S. W. Wells M. E. (1985). Reproductive Criteria of Beef Bulls during and after Exposure to Increased Ambient Temperature. J. Anim. Sci. 60 (2), 352–357. 10.2527/jas1985.602352x Mi H. Muruganujan A. Ebert D. Huang X. Thomas P. D. (2019). PANTHER Version 14: More Genomes, a New PANTHER GO-Slim and Improvements in Enrichment Analysis Tools. Nucleic Acids Res. 47 (D1), D419–D426. 10.1093/nar/gky1038 Mohammadi A. Alijani S. Rafat S. A. Abdollahi-Arpanahi R. (2020). Genome-Wide Association Study and Pathway Analysis for Female Fertility Traits in Iranian Holstein Cattle. Ann. Anim. Sci. 20 (3), 825–851. 10.2478/aoas-2020-0031 Mommersteeg M. T. M. Yeh M. L. Parnavelas J. G. Andrews W. D. (2015). Disrupted Slit-Robo Signalling Results in Membranous Ventricular Septum Defects and Bicuspid Aortic Valves. Cardiovasc. Res. 106 (1), 55–66. 10.1093/cvr/cvv040 Monteiro A. P. A. Tao S. Thompson I. M. T. Dahl G. E. (2016). In Utero heat Stress Decreases Calf Survival and Performance through the First Lactation. J. Dairy Sci. 99 (10), 8443–8450. 10.3168/jds.2016-11072 Mota L. F. M. Lopes F. B. Fernandes Júnior G. A. Rosa G. J. M. Magalhães A. F. B. Carvalheiro R. (2020). Genome-wide Scan Highlights the Role of Candidate Genes on Phenotypic Plasticity for Age at First Calving in Nellore Heifers. Sci. Rep. 10 (1), 6481. 10.1038/s41598-020-63516-4 Namous H. Peñagaricano F. Del Corvo M. Capra E. Thomas D. L. Stella A. (2018). Integrative Analysis of Methylomic and Transcriptomic Data in Fetal Sheep Muscle Tissues in Response to Maternal Diet during Pregnancy. BMC Genomics 19 (1), 123. 10.1186/s12864-018-4509-0 Ouellet V. Laporta J. Dahl G. E. (2020). Late Gestation Heat Stress in Dairy Cows: Effects on Dam and Daughter. Theriogenology 150, 471–479. 10.1016/j.theriogenology.2020.03.011 Peñagaricano F. Wang X. Rosa G. J. Radunz A. E. Khatib H. (2014). Maternal Nutrition Induces Gene Expression Changes in Fetal Muscle and Adipose Tissues in Sheep. BMC Genomics 15 (1), 1034. 10.1186/1471-2164-15-1034 Peterson D. G. Matitashvili E. A. Bauman D. E. (2003). Diet-Induced Milk Fat Depression in Dairy Cows Results in Increased Trans-10, Cis-12 CLA in Milk Fat and Coordinate Suppression of mRNA Abundance for Mammary Enzymes Involved in Milk Fat Synthesis. J. Nutr. 133 (10), 3098–3102. 10.1093/jn/133.10.3098 Pierce C. F. Speidel S. E. Coleman S. J. Enns R. M. Bailey D. W. Medrano J. F. (2020). Genome-wide Association Studies of Beef Cow Terrain-Use Traits Using Bayesian Multiple-SNP Regression. Livestock Sci. 232, 103900. 10.1016/j.livsci.2019.103900 Pimentel E. C. G. Bauersachs S. Tietze M. Simianer H. Tetens J. Thaller G. (2011). Exploration of Relationships between Production and Fertility Traits in Dairy Cattle via Association Studies of SNPs within Candidate Genes Derived by Expression Profiling. Anim. Genet. 42 (3), 251–262. 10.1111/j.1365-2052.2010.02148.x Ponsuksili S. Murani E. Schwerin M. Schellander K. Tesfaye D. Wimmers K. (2012). Gene Expression and DNA-Methylation of Bovine Pretransfer Endometrium Depending on its Receptivity after In Vitro-Produced Embryo Transfer. PLOS ONE 7 (8), e42402. 10.1371/journal.pone.0042402 Prasad S. S. Garg A. Agarwal A. K. (2011). Enzymatic Activities of the Human AGPAT Isoform 3 and Isoform 5: Localization of AGPAT5 to Mitochondria. J. Lipid Res. 52 (3), 451–462. 10.1194/jlr.M007575 Puig-Oliveras A. Ballester M. Corominas J. Revilla M. Estellé J. Fernández A. I. (2014). A Co-association Network Analysis of the Genetic Determination of Pig Conformation, Growth and Fatness. PLOS ONE 9 (12), e114862. 10.1371/journal.pone.0114862 Quinlan A. R. Hall I. M. (2010). BEDTools: A Flexible Suite of Utilities for Comparing Genomic Features. Bioinformatics 26 (6), 841–842. 10.1093/bioinformatics/btq033 R Core Team (2019). R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing. Available at: https://www.R-project.org/. Raeymaekers J. A. M. Chaturvedi A. Hablützel P. I. Verdonck I. Hellemans B. Maes G. E. (2017). Adaptive and Non-adaptive Divergence in a Common Landscape. Nat. Commun. 8 (1), 267. 10.1038/s41467-017-00256-6 Razzaghi A. Valizadeh R. Naserian A. A. Mesgaran M. D. Carpenter A. J. Ghaffari M. H. (2016). Effect of Dietary Sugar Concentration and sunflower Seed Supplementation on Lactation Performance, Ruminal Fermentation, Milk Fatty Acid Profile, and Blood Metabolites of Dairy Cows. J. Dairy Sci. 99 (5), 3539–3548. 10.3168/jds.2015-10565 Richards E. J. (2011). Natural Epigenetic Variation in Plant Species: a View from the Field. Curr. Opin. Plant Biol. 14 (2), 204–209. 10.1016/j.pbi.2011.03.009 Robinson J. J. Ashworth C. J. Rooke J. A. Mitchell L. M. McEvoy T. G. (2006). Nutrition and Fertility in Ruminant Livestock. Anim. Feed Sci. Tech. 126 (3), 259–276. 10.1016/j.anifeedsci.2005.08.006 Robinson M. D. McCarthy D. J. Smyth G. K. (2010). edgeR: A Bioconductor Package for Differential Expression Analysis of Digital Gene Expression Data. Bioinformatics 26 (1), 139–140. 10.1093/bioinformatics/btp616 Russo V. Fontanesi L. Davoli R. Nanni Costa L. Cagnazzo M. Buttazzoni L. (2002). Investigation of Candidate Genes for Meat Quality in Dry-Cured Ham Production: the Porcine Cathepsin B (CTSB ) and Cystatin B (CSTB ) Genes. Anim. Genet. 33 (2), 123–131. 10.1046/j.1365-2052.2002.00835.x Salleh M. S. Mazzoni G. Höglund J. K. Olijhoek D. W. Lund P. Løvendahl P. (2017). RNA-seq Transcriptomics and Pathway Analyses Reveal Potential Regulatory Genes and Molecular Mechanisms in High- and Low-Residual Feed Intake in Nordic Dairy Cattle. BMC Genomics 18 (1), 258. 10.1186/s12864-017-3622-9 Savolainen O. Lascoux M. Merilä J. (2013). Ecological Genomics of Local Adaptation. Nat. Rev. Genet. 14 (11), 807–820. 10.1038/nrg3522 Sejian V. Maurya V. P. Naqvi S. M. K. (2010). Adaptability and Growth of Malpura Ewes Subjected to thermal and Nutritional Stress. Trop. Anim. Health Prod. 42 (8), 1763–1770. 10.1007/s11250-010-9633-z Sevane N. Martínez R. Bruford M. W. (2018). Genome-wide Differential DNA Methylation in Tropically Adapted Creole Cattle and Their Iberian Ancestors. Anim. Genet. 50 (1), 15–26. 10.1111/age.12731 Shin D. Lee C. Park K.-D. Kim H. Cho K.-h. (2017). Genome-association Analysis of Korean Holstein Milk Traits Using Genomic Estimated Breeding Value. Asian-australas J. Anim. Sci. 30 (3), 309–319. 10.5713/ajas.15.0608 Skibiel A. L. Dado-Senn B. Fabris T. F. Dahl G. E. Laporta J. (2018). In Utero exposure to thermal Stress Has Long-Term Effects on Mammary Gland Microstructure and Function in Dairy Cattle. PLOS ONE 13 (10), e0206046. 10.1371/journal.pone.0206046 Skibiel A. L. Peñagaricano F. Amorín R. Ahmed B. M. Dahl G. E. Laporta J. (2018). In Utero Heat Stress Alters the Offspring Epigenome. Sci. Rep. 8 (1), 14609. 10.1038/s41598-018-32975-1 Storey J. D. Bass A. J. Dabney A. Robinson D. (2015). Qvalue: Q-Value Estimation for False Discovery Rate Control. R package version 2.8.0. . Http://Github.Com/Jdstorey/Qvalue. Sugimoto M. Fujikawa A. Womack J. E. Sugimoto Y. (2006). Evidence that Bovine Forebrain Embryonic Zinc finger-like Gene Influences Immune Response Associated with Mastitis Resistance. Proc. Natl. Acad. Sci. 103 (17), 6454–6459. 10.1073/pnas.0601015103 Suzuki M. M. Bird A. (2008). DNA Methylation Landscapes: Provocative Insights from Epigenomics. Nat. Rev. Genet. 9 (6), 465–476. 10.1038/nrg2341 Taberlet P. Valentini A. Rezaei H. R. Naderi S. Pompanon F. Negrini R. (2008). Are Cattle, Sheep, and Goats Endangered Species? Mol. Ecol. 17 (1), 275–284. 10.1111/j.1365-294X.2007.03475.x Takeuchi K. Reue K. (2009). Biochemistry, Physiology, and Genetics of GPAT, AGPAT, and Lipin Enzymes in Triglyceride Synthesis. Am. J. Physiology-Endocrinology Metab. 296 (6), E1195–E1209. 10.1152/ajpendo.90958.2008 Tao S. Bubolz J. W. do Amaral B. C. Thompson I. M. Hayen M. J. Johnson S. E. (2011). Effect of Heat Stress during the Dry Period on Mammary Gland Development. J. Dairy Sci. 94 (12), 5976–5986. 10.3168/jds.2011-4329 Vallée M. Beaudry D. Roberge C. Matte J. J. Blouin R. Palin M.-F. (2003). Isolation of Differentially Expressed Genes in Conceptuses and Endometrial Tissue of Sows in Early Gestation1. Biol. Reprod. 69 (5), 1697–1706. 10.1095/biolreprod.103.019307 Waggott D. Haider S. Boutros P. C. (2017). Bedr: Genomic Region Processing Using Tools Such as “BEDTools”, “BEDOPS” and “Tabix”. R package version 1.0.4. Https://CRAN.R-Project.Org/Package=bedr. Wang K. H. Brose K. Arnott D. Kidd T. Goodman C. S. Henzel W. (1999). Biochemical Purification of a Mammalian Slit Protein as a Positive Regulator of Sensory Axon Elongation and Branching. Cell 96 (6), 771–784. 10.1016/S0092-8674(00)80588-7 Wang X. Lan X. Radunz A. E. Khatib H. (2015). Maternal Nutrition during Pregnancy Is Associated with Differential Expression of Imprinted Genes and DNA Methyltranfereases in Muscle of Beef Cattle Offspring1. J. Anim. Sci. 93 (1), 35–40. 10.2527/jas.2014-8148 Weir B. S. Cockerham C. C. (1984). Estimating F -Statistics for the Analysis of Population Structure. Evolution 38 (6), 1358–1370. 10.2307/240864110.1111/j.1558-5646.1984.tb05657.x Xiang R. MacLeod I. M. Bolormaa S. Goddard M. E. (2017). Genome-wide Comparative Analyses of Correlated and Uncorrelated Phenotypes Identify Major Pleiotropic Variants in Dairy Cattle. Sci. Rep. 7 (1), 9248. 10.1038/s41598-017-09788-9 Yang J. Jiang J. Liu X. Wang H. Guo G. Zhang Q. (2016). Differential Expression of Genes in Milk of Dairy Cattle during Lactation. Anim. Genet. 47 (2), 174–180. 10.1111/age.12394 Yang J. Li W.-R. Lv F.-H. He S.-G. Tian S.-L. Peng W.-F. (2016). Whole-Genome Sequencing of Native Sheep Provides Insights into Rapid Adaptations to Extreme Environments. Mol. Biol. Evol. 33 (10), 2576–2592. 10.1093/molbev/msw129 Zeder M. A. (2008). Domestication and Early Agriculture in the Mediterranean Basin: Origins, Diffusion, and Impact. Proc. Natl. Acad. Sci. 105 (33), 11597–11604. 10.1073/pnas.0801317105 Zhang X. Zhang S. Ma L. Jiang E. Xu H. Chen R. (2017). Reduced Representation Bisulfite Sequencing (RRBS) of Dairy Goat Mammary Glands Reveals DNA Methylation Profiles of Integrated Genome-wide and Critical Milk-Related Genes. Oncotarget 8 (70), 115326–115344. 10.18632/oncotarget.23260 Zhao C. Gui L. Li Y. Plath M. Zan L. (2015). Associations between Allelic Polymorphism of the BMP Binding Endothelial Regulator and Phenotypic Variation of Cattle. Mol. Cell Probes 29 (6), 358–364. 10.1016/j.mcp.2015.09.007 Zhou D. a. Wei Z. Deng S. Wang T. Zai M. Wang H. (2013). SASH1 Regulates Melanocyte Transepithelial Migration through a Novel Gαs-SASH1-IQGAP1-E-Cadherin Dependent Pathway. Cell Signal. 25 (6), 1526–1538. 10.1016/j.cellsig.2012.12.025 Zhou X. Lindsay H. Robinson M. D. (2014). Robustly Detecting Differential Expression in RNA Sequencing Data Using Observation Weights. Nucleic Acids Res. 42 (11), e91. 10.1093/nar/gku310 Zoldoš V. Biruš I. Muratović E. Šatović Z. Vojta A. Robin O. (2018). Epigenetic Differentiation of Natural Populations of Lilium Bosniacum Associated with Contrasting Habitat Conditions. Genome Biol. Evol. 10 (1), 291–303. 10.1093/gbe/evy010