# Genome-Wide Analysis Reveals Selection Signatures Involved in Meat Traits and Local Adaptation in Semi-Feral Maremmana Cattle

Slim Ben-Jemaa, Gabriele Senczuk, Elena Ciani, Roberta Ciampolini, Gennaro Catillo, Mekki Boussaha, Fabio Pilla, Baldassare Portolano, Salvatore Mastrangelo

<figure><img src="image_1.png" /></figure>

## To cite this version:

Slim Ben-Jemaa, Gabriele Senczuk, Elena Ciani, Roberta Ciampolini, Gennaro Catillo, et al. Genome-Wide Analysis Reveals Selection Signatures Involved in Meat Traits and Local Adaptation in Semi-Feral Maremmana Cattle. Frontiers in Genetics, 2021, 10.3389/fgene.2021.675569. [hal-03210766](https://hal.inrae.fr/hal-03210766)

HAL Id: hal-03210766

https://hal.inrae.fr/hal-03210766

Submitted on 28 Apr 2021

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

Check for
updates

# Genome-Wide Analysis Reveals Selection Signatures Involved in Meat Traits and Local Adaptation in Semi-Feral Maremmana Cattle

Slim Ben-Jemaa¹, Gabriele Senczuk², Elena Ciani³, Roberta Ciampolini⁴, Gennaro Catillo⁵,
Mekki Boussaha⁶, Fabio Pilla², Baldassare Portolano⁷ and Salvatore Mastrangelo⁷*

¹Laboratoire des Productions Animales et Fourragères, Institut National de la Recherche Agronomique de Tunisie, University of Carthage, Ariana, Tunisia, ²Dipartimento di Agricoltura, Ambiente e Alimenti, University of Molise, Campobasso, Italy, ³Dipartimento di Bioscienze, Biotecnologie e Biofarmaceutica, University of Bari "Aldo Moro", Bari, Italy, ⁴Dipartimento di Scienze Veterinarie, University of Pisa, Pisa, Italy, ⁵Consiglio per la Ricerca in Agricoltura e l'Analisi dell'Economia Agraria (CREA), Centro di Ricerca Zootecnia e Acquacoltura, Lodi, Italy, ⁶INRAE, AgroParisTech, University of Paris Saclay, Saint Aubin, France, ⁷Dipartimento di Scienze Agrarie, Alimentari e Forestali, University of Palermo, Palermo, Italy

OPEN ACCESS

**Edited by:**
Emiliano Lasagna,
University of Perugia, Italy

**Reviewed by:**
Francisco Javier Navas González,
University of Cordoba, Spain
Giustino Gaspa,
University of Turin, Italy
Nuno Carolino,
Instituto Nacional Investigacao
Agraria e Veterinaria (INIAV), Portugal

*Correspondence:
Salvatore Mastrangelo
salvatore.mastrangelo@unipa.it

**Specialty section:**
This article was submitted to Livestock Genomics, a section of the journal Frontiers in Genetics

Received: 03 March 2021

Accepted: 07 April 2021

Published: 28 April 2021

**Citation:**
Ben-Jemaa S, Senczuk G, Ciani E, Ciampolini R, Catillo G, Boussaha M, Pilla F, Portolano B and Mastrangelo S (2021) Genome-Wide Analysis Reveals Selection Signatures Involved in Meat Traits and Local Adaptation in Semi-Feral Maremmana Cattle. *Front. Genet.* 12:675569.
doi: 10.3389/fgene.2021.675569

The Maremmana cattle is an ancient Podolian-derived Italian breed raised in semi-wild conditions with distinctive morphological and adaptive traits. The aim of this study was to detect potential selection signatures in Maremmana using medium-density single nucleotide polymorphism array. Putative selection signatures were investigated combining three statistical approaches designed to quantify the excess of haplotype homozygosity either within (integrated haplotype score, iHS) or among pairs of populations (Rsb and XP-EHH), and contrasting the Maremmana with a single reference population composed of a pool of seven Podolian-derived Italian breeds. Overall, the three haplotype-based analyses revealed selection signatures distributed over 19 genomic regions. Of these, six relevant candidate regions were identified by at least two approaches. We found genomic signatures of selective sweeps spanning genes related to mitochondrial function, muscle development, growth, and meat traits (SCIN, THSD7A, ETV1, UCHL1, and MYOD1), which reflects the different breeding schemes between Maremmana (semi-wild conditions) and the other Podolian-derived Italian breeds (semi-extensive). We also identified several genes linked to Maremmana adaptation to the environment of the western-central part of Italy, known to be hyperendemic for malaria and other tick-borne diseases. These include several chemokine (C-C motif) ligand genes crucially involved in both innate and adaptive immune responses to intracellular parasite infections and other genes playing key roles in pulmonary disease (HEATR9, MMP28, and ASIC2) or strongly associated with malaria resistance/susceptibility (AP2B1). Our results provide a glimpse into diverse selection signatures in Maremmana cattle and can be used to enhance our understanding of the genomic basis of environmental adaptation in cattle.

**Keywords:** local cattle breeds, selection signatures, diversity, candidate genes, environmental adaptation
INTRODUCTION

In livestock species, domestication followed by breed formation and the subsequent artificial selection for economic and morphological traits have shaped genomic variation, driving the formation of detectable signatures on the genome of these populations (e.g., Stella et al., 2010; Qanbari et al., 2011; Rothammer et al., 2013; Signer-Hasler et al., 2017; Mastrangelo et al., 2020a). With more than 30 officially recognized local cattle breeds, Italy can be considered as one of the most important hotspots of cattle diversity (Mastrangelo et al., 2018). Several Italian breeds such as Marchigiana, Romagnola, Chianina, Calvana, Podolica, Cinisara, Modicana, and Maremmana belong to the so-called “Podolian” group (hereinafter referred to as “Podolian-derived Italian breeds”; Mastrangelo et al., 2020b). This group includes several European breeds with common phenotypic traits such as a gray coat color and long horns (Di Lorenzo et al., 2018) considered to be ancestral. Indeed, their evolutionary history seems to be linked to human Neolithic and post-Neolithic migration routes (Senczuk et al., 2020).

The Maremmana is an ancient cattle breed that, unlike the other Podolian-derived Italian breeds, is raised in semi-feral conditions (Figure 1), whose origin can be traced back to the Etruscan era when local would have been crossed with Podolian cattle (Maretto et al., 2012). This breed has been reared for hundreds of years in the lowlands and on the hilly areas of Maremma, a territory bordering the Tyrrhenian Sea and including much of south-western Tuscany and part of northern Lazio, in western central Italy. This region is characterized by an arid climate and poor quality pastures (Food and Agriculture Organization, 1991). The Maremma territory is also known to have been hyperendemic for malaria since the 60s of the last century (Romi et al., 2012). Until the 20th century, Maremmana

cattle ran semi-wild driven by mounted herdsmen called “Vaccari” or “Butteri” (Ritchie, 2009). This breed is well-known to be resistant to several diseases such as tuberculosis, parasitic diseases (e.g., piroplasmosis), and foot and skin diseases (Food and Agriculture Organization, 1991). The Maremmana cattle, compared with other Podolian-derived Italian breeds, still display many ancestral features such sexual dimorphism, good maternal qualities, disease resistance, and very long and distinctively lyre-shaped horns. The breed has also a robust appearance characterized by a massive skeletal structure, solid legs, and hard hooves, longevity and adaptation to harsh conditions and poor-quality feed. Because of mechanization in agriculture and drainage of the marshy area, the population size of Maremmana has fallen from ~288,000, in 1944 (Consiglio Nazionale delle Ricerche, 1983) to roughly 44,000 in 1983. Since 70s, the breed was used as “maternal line” in crossbreeding programs with local and cosmopolite beef breeds. The indiscriminate and unregulated adoption of crossbreeding induced a constant and progressive depletion of pure animals to few 1,000 females (Fioretti et al., 2020). By the end of the 1980’s, this breed was nearly extinct (Food and Agriculture Organization, 2002). Later on, following the implementation of agri-environmental measures and the introduction of incentives to protect endangered livestock populations, the trend started to reverse and the number of individuals belonging to Maremmana population increased (Catillo et al., 2018). In fact, about 11,800 animals are recorded in the herdbook in 2019, including 6,853 cows and 213 bulls reared in 250 herds for beef production.<fnref n="1" />

Previous studies have investigated genetic diversity and population structure of Maremmana and its relationship with other breeds (e.g., Ajmone-Marsan et al., 2002; Moioli et al., 2004; Pariset et al., 2010; Maretto et al., 2012; Catillo et al., 2018; Biscarini et al., 2020; Senczuk et al., 2020), but to the best of our knowledge, there is no information on trait-specific selection signatures for this breed. Such information may contribute to improve our understanding of how selection affects livestock phenotypes, thus allowing important insights into the evolutionary history and the genetic architecture of the selected traits.

Several studies aiming at identifying genomic regions putatively under selection in cattle have been conducted by contrasting breeds with different production aptitudes (dairy and meat; e.g., Stella et al., 2010; Bomba et al., 2015; Zhao et al., 2015). Other studies have shown that it is possible to identify selection signatures by comparing breeds with similar selection goals (Sorbolini et al., 2015). The rationale behind this approach is that even though such breeds share selection goals, every breed has, in the course of time, been subjected to specific selection pressures, which might affect different set of genes.

The aim of this study was to map selection signatures in Maremmana cattle that may be linked to peculiar phenotypic characteristics. For this purpose, three statistical approaches designed to quantify the excess of haplotype homozygosity (relative to neutral expectations) were applied: (a) integrated haplotype score (iHS); (b) Rsb [standardized log-ratio of the integrated site-specific Extended Haplotype Homozygosity between pairs of populations]; and (c) XP-EHH (cross-population EHH).

<figure><img src="image_4.png" /><figcaption>FIGURE 1 | A typical herd of Maremmana cattle.</figcaption></figure>

<footnote n="1">¹http://www.anabic.it/index1.htm</footnote>
For Rsb and *XP-EHH* computation, the excess of long haplotypes quantified within Maremmana was contrasted with that estimated from seven other Podolian-derived Italian breeds (considered as the meta-population). *iHS* (Voight et al., 2006) detects partial selective sweeps, while Rsb (Tang et al., 2007) and *XP-EHH* (Sabeti et al., 2007) detect selected alleles that have risen (in frequency) to near fixation.

# MATERIALS AND METHODS

## Genotyping Data, Quality Control, and Relatedness Filtering

Genotypes at loci in the Bovine SNP50K v2 BeadChip (Illumina Inc., San Diego, CA, United States) for 146 Maremmana (135 females and 11 males) and 174 pooled samples from seven other Podolian-derived Italian breeds (Marchigiana, Romagnola, Chianina, Calvana, Podolica, Cinisara and Modicana) were retrieved from previously published studies (Mastrangelo et al., 2018; Biscarini et al., 2020; Senczuk et al., 2020). In general, methods for the detection of selection signatures are designed to analyze non-related animals (Zhao et al., 2015). For this purpose, a relatedness test between individuals within each population was first performed using PLINK (Purcell et al., 2007). The software calculates a variable called PIHAT reflecting extended haplotypes shared between distantly related individuals. One individual from any pairs showing a PIHAT score ≥0.15 was removed from further analysis. After relatedness filtering, 71 Maremmana and 163 Italian Podolian-derived individuals were retained for analyses. Further quality control (QC) filtering steps to minimize the number of false-positive calls were applied. These filters included the removal of markers with minor allele frequency (MAF) ≤ 0.01 and a call rate ≤ 0.95. A total of 34,973 single nucleotide polymorphisms (SNPs) were kept after these QC steps. Moreover, samples with more than 5% of missing genotypes were also discarded from further analyses.

Because of uncertainty in the identification of the ancestral allelic state (see the section “Identification of selection signatures” below) of 3,223 SNPs in *iHS* computation, we applied this statistic to a second dataset with 31,750 SNPs. The two datasets were then subjected to linkage disequilibrium (LD) pruning using the default parameters of PLINK (SNP window size: 50, step 5 SNPs, r<fnref n="2" />=0.5) resulting in the removal of 9,517 and 3,265 SNPs from the *iHS* and the *Rsb/XP-EHH* datasets, respectively. Overall, 22,233, 31,708 and 31,708 markers were used for *iHS*, *Rsb*, and *XP-EHH* computations, respectively.

## Genetic Relationships

To investigate the existence of population substructure for the breeds of the study, identity-by-state (IBS) alleles shared by all pairs of individuals were calculated using PLINK (Purcell et al., 2007) and graphically represented by multidimensional scaling (MDS).

Moreover, we computed a matrix of distances between all pairs of individuals using the ape R package (Paradis et al., 2004). The genetic distance between two individuals was defined as the number of SNPs for which they differ. A neighbor-joining

tree was then computed based on the resulting distance matrix using the phyclust R package (Chen, 2011). Both analyses were carried out using the dataset containing 31,708 SNPs and 234 individuals (71 Maremmana and 163 Italian Podolian-derived).

## Identification of Selection Signatures

The three EHH-based metrics were computed using the rehh package in R (Gautier and Vitalis, 2012). Haplotypes were reconstructed from the genotyped SNPs using fastPHASE (Scheet and Stephens, 2006). For each chromosome, the following options were used: -T10 -Ku60 -Kl10 -Ki10, where: -T is the number of random starts of the expectation-maximization (EM) algorithm used by fastPHASE; and -Ku is the upper limit for the number of considered clusters; -Kl is the lower limit for the number of considered clusters; -Ki is the interval between values for number of clusters.

The *iHS* statistic measures the amount of EHH for a given SNP along the ancestral allele relative to the derived allele. In this study, for each SNP, the ancestral allele was inferred as the most common allele within three outgroup species including yak, buffalo, and sheep. For each SNP, an *iHS* score was computed following Voight et al. (2006). SNP scores were further transformed into two-sided *p*-values: $piHS = -\log10[1-2|\Phi(iHS)-0.5|]$, where $\Phi(x)$ represents the Gaussian cumulative distribution function (assuming that the *iHS* values are normally distributed, which is expected under neutrality). We used $-\log10(p-value) = 2.5$ as a threshold to define significant *iHS* values. Because it has been demonstrated that it is more powerful to look for windows of consecutive SNPs that contain numerous extreme *iHS* scores (Voight et al., 2006), *iHS* candidate regions were defined as those containing at least three neighboring SNPs exceeding the aforementioned threshold within 2-Mb sliding windows, overlapping by 20 Kb.

The Rsb and *XP-EHH* scores were computed between Maremmana and the other Podolian-derived breeds following Sabeti et al. (2007) and Tang et al. (2007), respectively. As in *iHS*, assuming that Rsb and *XP-EHH* values are normally distributed, SNP scores were further transformed into two-sided *p*-values: $pRsb = -\log10[1-2|\Phi(Rsb)-0.5|]$ and $pXP-EHH = -\log10[1-2|\Phi(XP-EHH)-0.5|]$. Rsb and *XP-EHH* candidate regions were defined as those containing at least four neighboring SNPs exceeding the $-\log10(p-value) = 2.5$ threshold within 1-Mb sliding windows, overlapping by 10 Kb.

## Functional Characterization of Regions Identified as Under Selection

Genomic regions detected by at least two EHH-based methods were interrogated for genes annotated to the *Bos taurus* genome assembly ASR-UCD1.2 using Genome Data Viewer provided by NCBI. Database for Annotation, Visualization and Integrated Discovery (DAVID) software was used for functional enrichment analysis of the list of genes retrieved.<fnref n="2" /> We used an adjusted Benjamini-corrected *p*-value = 0.05 as a threshold to define significantly enriched Gene Ontology (GO) biological and

<footnote n="2">²https://david.ncifcrf.gov/</footnote>
functional processes. Moreover, to investigate the biological function and the phenotypes that are known to be affected by each annotated gene, we conducted a comprehensive literature search, including information from other species.

## RESULTS AND DISCUSSION

The present study aimed at the detection of genomic regions that had been differentially selected between Maremmana and other Podolian-derived Italian breeds, which can provide clues to potential new targets of natural and artificial selection. For this purpose, a genomic scan for selective sweeps was performed by contrasting the extended haplotype homozygosity profiles between Maremmana and a single reference population composed of a pool of seven Podolian-derived Italian breeds (Rsb and XP-EHH). The third complementary EHH-based method, iHS, was applied because it is known to be sensitive to ongoing within-population recent selection (Voight et al., 2006).

### Genetic Relationships

Multidimensional scaling analysis on IBS allele-sharing values revealed a clear separation between Maremmana and the other

Podolian-derived Italian breeds along the first dimension (C1; **Figure 2**). This clustering trend was confirmed by the neighbor joining (NJ) tree of individuals, which revealed two major clades that coincide with Maremmana and the rest of the Podolian-derived populations (**Supplementary Figure S1**), in agreement with previous studies (Senczuk et al., 2020; Zsolnai et al., 2020). It is also noteworthy that, both MDS and NJ tree revealed additional substructure of Maremmana suggesting the presence of individuals from several sub-populations with different allele frequencies. This genetic background is likely the result of bottleneck in the population, probably promoted by the overuse of few influential sires as parents of the next generation in the early 1990s (Porter et al., 2016), or the consequence of feral breeding. Moreover, for Maremmana breed, natural mating system is the common practice and the exchange of bulls among herds is quite unusual. This leads to an increase of inbreeding within sub-populations.

### Selection Signatures Detected by the Within-Population Test

The use of *iHS* test aims to identify recent and incomplete selective sweeps (Voight et al., 2006). This approach revealed a total of 103 SNPs putatively under selection in Maremmana breed.

<figure><img src="image_5.png" /><figcaption>**FIGURE 2** | Genetic relationship defined with multidimensional scaling (MDS) analysis between Maremmana and Podolian-derived breeds.</figcaption></figure>
These markers defined a total of three regions, two located on BTA04 (at positions 18.44-22.80 and 24.62-28.60 Mb) and one on BTA12 (44.98-47.12 Mb; **Figure 3**), spanning a total of 28 and 25 known and uncharacterized genes (LOC), respectively (**Supplementary Table S1**). The region with the strongest signal (eight SNPs exceeding the significance threshold) is located on the BTA04 at position 18.44-22.80 Mb and spans six known genes (*THSD7A*, *TMEM106B*, *VWDE*, *SCIN*, *ARL4A*, and *ETV1*). *SCIN* and *THSD7A* genes are involved in several biological processes (BP) related to actin, an important contributor to the contractile property of muscle, such as filament capping (GO: 0051693), and severing (GO: 0051014), polymerization or depolymerization (GO:0008154) and cytoskeleton reorganization (GO:0031532). Actin is the building block of thin filaments of skeletal muscle (Lonergan et al., 2010), and together with myosin, it can affect meat quality post-mortem. A third candidate for genes under selection in the region on BTA04 is *ETV1*, a member of the PEA3 group of transcription factors, with functions related, inter alia, to muscle organ development (GO:0007517). *ETV1* is one of the

strongest candidate genes for marbling score in Korean cattle (Ryu and Lee, 2016). Likewise, in a previous study analyzing the muscle transcriptome in pigs, *ETV1* was found among the top 20 differentially expressed genes between Duroc and Landrace breeds in a feed efficiency context (Carmelo and Kadarmideen, 2020).

## Selection Signatures Detected by Between-Populations Tests

The clear separation between Maremmana and the other Podolian-derived Italian breeds, in conjunction with Maremmana membership to the Podolian group suggest that between-population tests offer a powerful approach to unravel potential selection differences between them. The choice of contrasting the EHH patterns of the same haplotype between Maremmana and the seven other Podolian-derived Italian breeds relied on the distinctive primitive features of the former and its long-term adaptability to the specific environmental conditions of the Maremma plain, which houses several protozoal parasites affecting cattle (Toma et al., 2015).

<figure><img src="image_6.png" /><figcaption>**FIGURE 3** | Manhattan plot of the genome-wide integrated haplotype score (iHS) analysis for Maremmana breed. Horizontal dashed line mark the significance threshold applied to detect the outlier single nucleotide polymorphisms (SNPs) –log10 (*p*-value) = 2.5.</figcaption></figure>
This region was also known to be hyper-endemic for malaria until the mid-20s (Boccolini et al., 2012).

We applied *Rsb* and *XP-EHH* tests to detect potential selective sweeps that are fixed (or nearly fixed) in Maremmana but remain segregating in the other Podolian-derived Italian breeds or vice versa.

In total, 396 and 421 outlier SNPs were revealed by *Rsb* and *XP-EHH* tests, respectively (Figures 4A,B), with the two approaches defining eight candidate regions each (*Supplementary Tables* S2 and S3) and spanning 141 and 219 known protein coding genes. Among these, four genomic regions located on BTA06, BTA15, BTA19, and BTA24 were identified by both approaches (*Rsb* and *XP-EHH*; Table 1). The outlier window on BTA06, spanning 5.54 Mb (at position 58.58–64.12 Mb), showed the strongest signal with 46 and 85% of the SNPs exceeding the significance threshold for the *Rsb* and the *XP-EHH* tests, respectively (Figure 4). The distinctive peak on BTA06 partially overlaps with a run of homozygosity (ROH) island previously detected in Maremmana (Biscarini et al., 2020)

and spanning four genes (*KCTD8*, *YIPF7*, *GUF1*, and *GNPDA2*) associated with carcass traits in a composite beef cattle breed (Hay and Roberts, 2018). The *GNPDA2* gene has a potential role in the regulation of body weight, fat, and energy metabolism in chickens (Ouyang et al., 2016) and was associated with adipose tissue accumulation and obesity in humans (Renström et al., 2009). Therefore, this region potentially contributes to the distinctive characteristics of Maremmana related to growth, carcass, and meat traits. The region on BTA06 also showed nine genes (*CHRNA9*, *LIAS*, *N4BP2*, *PDS5A*, *RBM47*, *RHOH*, *SMIM14*, *UBE2K*, and *UGDH*) that overlaps with a ROH island in fast-growing meat commercial hybrid Turkey (Strillacci et al., 2020). Notably, *CHRNA9* gene is also associated with meat quality traits in pigs (Velez-Irizarry et al., 2019). Moreover, the selection signature on BTA06 showed several genes related to mitochondrial functions. For instance, the protein encoded by *LIAS* gene localizes in mitochondrion and is involved in response to oxidative stress (GO:0006979). *GUF1* encodes a GTPase that triggers back-translocation of the elongating ribosome during

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /><figcaption>**FIGURE 4** | Manhattan plots showing the results of (A) *Rsb* and (B) cross-population extended haplotype homozygosity (*XP-EHH*) test analyses in the comparison between Maremmana and Italian Podolian-derived breeds. Horizontal dashed lines mark the significance threshold applied to detect the outlier SNPs –log10 (*p*-value) = 2.5.</figcaption></figure>
<table><caption>TABLE 1 | Overlapped genomic regions identified in at least two approaches in Maremmana cattle breed.</caption><thead><tr><th>BTA</th><th>Start (Mb)</th><th>End (Mb)</th><th>Approaches</th><th>Genes</th></tr></thead><tbody><tr><td>4</td><td>19.31</td><td>22.10</td><td>iHS, Rsb and XP-EHH</td><td>THSD7A, LOC101902990, TMEM106B, LOC112446322, LOC101906705, LOC104971960, SCIN, LOC100297402, ARL4A, LOC112446323, LOC781728, LOC112446324, LOC112446494, ETV1</td></tr><tr><td>6</td><td>58.58</td><td>64.12</td><td>Rsb and XP-EHH</td><td>LIAS, UGDH, SMIM14, UBE2K, LOC112447195, PDS5A, TRNAG-CCC, LOC112447155, N4BP2, RHOH, CHRNA9, LOC112447071, RBM47, LOC101901948, LOC104968867, NSUN7, APBB2, LOC613534, UCHL1, LIMCH1, LOC112447073, PHOX2B, TMEM33, LOC784473, SLC30A9, LOC782858, BEND4, LOC101903036, LOC107132567, LOC112447075, SHISA3, ATP8A1, LOC527955, GRXCR1, LOC112447156, LOC112447157, LOC101906152, KCTD8, YIPF7, GUF1, GNPDA2, LOC104968873, LOC112447181, LOC112447015</td></tr><tr><td>12</td><td>46.31</td><td>47.12</td><td>iHS and Rsb</td><td>DACH1</td></tr><tr><td>15</td><td>33.01</td><td>34.80</td><td>Rsb and XP-EHH</td><td>LOC107131357, UBASH3B, LOC107133160, CRTAM, JHY, LOC107133161, BSX, LOC107133162, HSPA8, LOC112441683, LOC112441682, CLMP, GRAMD1B, MIR2313, SCN3B, ZNF202, LOC112441623, LOC100847729, SAAL1, TPH1, SERGEF, KCNC1, MYOD1</td></tr><tr><td>19</td><td>14.02</td><td>17.36</td><td>Rsb and XP-EHH</td><td>HNF1B, LOC101906072, LOC104969035, LOC101906153, HEATR6, LOC100847618, LOC515676, LOC525415, LOC100462699, LOC100848100, WFDC18, LOC100296618, LOC100847724, CCL4, LOC107131498, CCL3, LOC616364, LOC100297044, LOC112442607, LOC504773, LOC508666, CCL14, CCL16, CCL5, HEATR9, TAF15, MMP28, C19H17orf50, GAS2L2, RASL10B, AP2B1, PEX12, LOC112442809, SLFN14, LOC112442608, SLFN11, LOC101907813, LOC100848263, LOC112442609, LOC112442610, LOC112442826, LOC112442590, LOC112442786, UNC45B, NLE1, FNDC8, RAD51D, MIR2331, RFFL, LIG3, LOC104969067, ZNF830, CCT6B, LOC112442567, LOC112442798, TMEM132E, LOC789304, CCL1, CCL8, CCL11, CCL2, LOC112442773, ASIC2, LOC104969069, LOC112442874, SPACA3, TMEM98, MYO1D</td></tr><tr><td>24</td><td>39.70</td><td>41.54</td><td>Rsb and XP-EHH</td><td>ARHGAP28, LAMA1, LOC101904673, LRRC30, LOC781276, PTPRM, MIR22850-4, RAB12, MTCL1, LOC107131785, LOC100849069, NDUFV2, ANKRD12</td></tr></tbody></table>

mitochondrial protein synthesis,<fnref n="3" /> while *UCHL1* is involved in regulation of mitochondrial content and function. Mouse model with skeletal muscle specific knockout of *UCHL1* caused a significant reduction in key proteins that are involved in mitochondrial oxidative phosphorylation in soleus muscles (Gao et al., 2020). This gene is also associated with several other processes such as adult walking behavior (GO:0007628) and eating behavior (GO:0042755). This distinctiveness between Maremmana and the other Podolian-derived Italian breeds is

believed to be due to differences in mitochondria-associated gene pathways. Indeed, meat quality of beef is likely to be impacted by breed differences in postmortem mitochondrial respiratory activity within the muscle tissue (Ramos et al., 2020). Furthermore, previous studies have suggested a link between muscle mitochondrial function and feed efficiency in several species, such as chickens (Ojano-Dirain et al., 2005), beef cattle (Fernandez et al., 2020), and pigs (Carmelo and Kadarmideen, 2020). Overall, several of our candidate genes affect multiple, apparently unrelated, phenotypes (pleiotropic effect), which represents one of the difficulties in identifying phenotypic targets of selection.

<footnote n="3">³https://www.genecards.org/cgi-bin/carddisp.pl?gene=GUF1</footnote>
# Relevant Selection Signatures of Distinctiveness in Maremmana Cattle Compared to Podolian-Derived Italian Breeds

A common constraint of selection signature detection methods is the identification of false signals that do not reflect selection (Fijarczyk and Babik, 2015); rather they arise from heterogeneous recombination rates along the genome or owing to several demographic processes such as migration, expansions, and bottlenecks (Vitti et al., 2013). Using multiple approaches and intersecting signals between at least two of them is recommended to reduce the number of false positives (de Villemereuil et al., 2014; Mastrangelo et al., 2020b). Of particular note is that, unlike previous studies (e.g., Bahbahani et al., 2015), we found a high congruence among the three EHH-based methods. This is particularly true for the *iHS* statistic, where two among the three outlier windows identified by this test, were also detected by the *Rsb* statistic (Table 1). To identify additional reliable candidate genes under differential selection between Maremmana and the other Podolian-derived Italian breeds, we focused on the six regions that were detected by at least two methods (Table 1).

## Genes Associated With Meat Quality Traits

In addition to the aforementioned genes in the selection signature on BTA04 (*ETV1*, *SCIN*, and *THSD7A*), we identified other candidate genes involved in several BP related to muscle compartments, and many of them were previously shown to be associated with meat quality traits. For instance, *MYOD1* on BTA15 (also known as the myogenic differentiation one gene), plays a key role in skeletal muscle cell differentiation (GO:0035914). Polymorphisms within this gene are associated with carcass and meat quality traits in beef cattle (Du et al., 2013). Furthermore, this gene overlapped with a selection region for intramuscular fat and backfat thickness in pig (Kim et al., 2015; Li et al., 2016). Likewise, the *LIMCH1* gene on BTA06 encodes an actin stress fibers-associated protein that activates non-muscle myosin IIA (Lin et al., 2017) and was shown to be under selection in Brahman cattle (Chen et al., 2020). Likewise, two other candidate genes on BTA19, *GAS2L2*, and *unc-45 myosin chaperone B* (*UNC45B*), were among possible candidate genes for 3 day *postmortem* collagen solubility in three Canadian beef cattle breeds (Lei, 2019). In particular, *GAS2L2* is mainly expressed in skeletal muscle. The protein encoded by this gene appears to cross-link microtubules and microfilaments in muscle (Stroud et al., 2014). *UNC45B* is essential for sarcomeric organization and muscle function in very distinct taxonomic groups spanning from *Caenorhabditis elegans* to humans (Donkervoort et al., 2020) and is among a group of genes that showed statistically significant differences in their expression between beef and dairy cattle (Ciecierska et al., 2020). Therefore, our findings partially reflect the different breeding schemes for meat quality traits between Maremmana and the other Podolian-derived Italian breeds. In fact, a previous study using high throughput mRNA sequencing revealed breed specific variations in gene expression of muscle tissues between Maremmana and Chianina

(Bongiorni et al., 2016), including two of our relevant candidate genes on BTA06, *UCHL1*, and *SHISA3*.

## Genes Associated With Immune Response and Local Adaptation

During hundreds of years, the Maremmana population has been exposed to vectors of several pathogens such as *Ixodes ricinus*, *Rhipicephalus bursa*, *Rhipicephalus turanicus*, and *Rhipicephalus sanguineus* (Toma et al., 2015) carrying hemoparasites such as Babesia and Theileria (Toma et al., 2017). Constant exposure to these parasites is expected to exert a selective pressure on Maremmana immunity-related genes. Accordingly, when we conducted enrichment analysis of the 78 protein coding genes identified in at least two approaches (Table 1), we found significant (Benjamini-corrected *p*-value < 0.05) functional enrichment of several BP and molecular function (MF) terms related to a variety of immune responses (Supplementary Table S4). Lymphocyte chemotaxis (GO:0048247, *n* = 7, Benjamini-corrected *p*-value = 1.32 × 10⁻⁷), monocyte chemotaxis (GO:0002548, *n* = 7, Benjamini-corrected *p*-value = 4.09 × 10⁻⁷), cellular response to interferon-gamma (GO:0071346, *n* = 7, Benjamini-corrected *p*-value = 1.09 × 10⁻⁶) and positive regulation of inflammatory response (GO:0050729, *n* = 7, Benjamini-corrected *p*-value = 2.16 × 10⁻⁶), were the most enriched BP terms. CCR chemokine receptor binding (GO:0048020, *n* = 7, Benjamini-corrected *p*-value = 3.26 × 10⁻⁸) and chemokine activity (GO:0008009, *n* = 7, Benjamini-corrected *p*-value = 4.07 × 10⁻⁷) were the most enriched terms under MF (Supplementary Table S4). Over-representation of GO functions playing a central role in the resistance to pathogens is due to the selection signal observed on the BTA19 harboring several chemokine (C-C motif) ligand genes (*CCL4*, *CCL3*, *CCL14*, *CCL16*, *CCL5*, *CCL1*, *CCL8*, *CCL11*, and *CCL2*). Chemokine ligands are crucially involved in both innate and adaptive responses, and help initiating the humoral immune defense (Shin et al., 2015). Through interaction with their receptors, they play a key role in attracting a diverse set of effector leukocytes to inflammatory sites and in recruiting neutrophils, monocytes/macrophages, dendritic cells (DC), and natural killer (NK) cells (Esche et al., 2005). These findings suggest that chemokine (C-C motif) ligand genes were differentially selected in the compared breeds. It seems that significant differences in cytokine production are behind the marked difference in resistance to challenge with *Theileria annulata* between *Bos taurus* and *Bos indicus* cattle subspecies (Glass et al., 2005). Likewise, by inspecting the relevant candidate regions that overlap between statistical methods, we also noted the presence of many genes with immune-related functions (*CRTAM*, *SLFN11*, *RAB12*, and *AP2B1*). The *CRTAM* gene on BTA15 is involved in adaptive immune response (GO:0002250), in positive regulation of cytokine secretion (GO:0001819) and positively regulates natural killer cell mediated cytotoxicity (GO:0045954). The *SLFN11* gene (on BTA19) belongs to the Schlafen protein family, which plays a prominent role in the induction of immune responses (Geserick et al., 2004) and in the regulation of viral replication (Li et al., 2012), while *RAB12* (on BTA24) is directly regulated in response to immune stimuli
(Efergan et al., 2016). The *AP2B1* (on BTA19) encodes a protein, which is a part of the AP2 coat assembly protein complex and links clathrin to receptors in the coated vesicles (Online Mendelian Inheritance in Man, 2003). Importantly, a SNP located in the human ortholog of this gene showed a high evidence for association with malaria susceptibility (Band et al., 2019). This is an intriguing outcome and is fairly compelling when considering a possible crucial role in the adaptation of Maremmana cattle to its malarial environment. Several of tick-borne diseases in the Maremma region (such as theileriosis and babesiosis) are also known to result into hypoxia, increased respiratory and heart rates (Razavi et al., 2015; Bal et al., 2016). Thus, it is not surprising that among our candidates, three genes located on the BTA19 (*HEATR9*, *MMP28*, and *ASIC2*), play key roles in pulmonary diseases. The *HEATR9* gene is highly upregulated in lung cells after influenza or respiratory syncytial virus (RSV) infections in humans (Stairiker et al., 2020). The *MMP28* is expressed by macrophages and regulates their recruitment to the lung, thus playing a role in pulmonary disease (Gharib et al., 2014), while *ASIC2* is expressed in the pulmonary arterial smooth muscle cells (Jernigan et al., 2009) and plays an important role in the regulation of pulmonary vascular reactivity and was shown to be under differential selection between African and Northern European Chickens (Fleming et al., 2017). Furthermore, two other candidate genes were also associated with adaptation to high altitude environmental condition (such as hypoxia) in several species: *SERGEF* on BTA15 (Li et al., 2013; Ma et al., 2019) and *NDUFV2* on BTA24 (Ning et al., 2010). *SERGEF* is also among the 48 differentially expressed genes from two breeds of chicken showing extreme differences in hypoxic adaptation reared in Tibet (Hu et al., 2020), whereas *NDUFV2* has been identified as a candidate gene possibly involved in pulmonary hypertension in Black Angus cattle (Newman et al., 2011). The *NDUFV2* gene has been also detected under selection pressure for lipid metabolism in Yaroslavl cattle, and related with local adaptation to survive under poor nutrition conditions in winter (Zinovieva et al., 2020). Additional genes located in the candidate regions on BTA15 (*HSPA8*, *SCN3B*) and BTA19 (*LAMA1*) point toward other aspects of environmental adaptation of Maremmana: *HSPA8* plays a role in thermoregulatory protective mechanisms in cattle and buffalo under tropical environments (Kumar et al., 2015), and associated with cellular response to heat stress in goats (Mohanarao et al., 2014); *SCN3B* is involved in positive regulation of heart rate (GO:0010460) and sensory perception of pain (GO:0019233). This gene is among the genes differentially expressed in response to hyperthermia induced apoptosis in broiler chicken (Coble et al., 2014). Finally, *LAMA1* has been identified as one of candidate genes implicated for environmental adaptation in indigenous Chinese cattle (Gao et al., 2017).

As expected, these findings support a model of complex integrated genomic circuits involving multiple pathways characterizing response mechanisms to environmental challenges in Maremmana. Yet, it seems unlikely that all the genes cited above are a target of selection. It is likely that only some of these genes may be associated with environmental adaptation in Maremmana, while the others would be indirectly influenced by natural selection through genetic hitchhiking.

## Genes Associated With Morphological Traits

Between population tests are expected to identify candidate genomic regions harboring genes related to some distinctive morphological primitive characteristics of Maremmana (e.g., horn shape and size, sexual dimorphism). However, our candidate regions harbored few morphology-related genes. *Notchless homolog 1* (*NLE1*; BTA19) is the only candidate gene directly involved in skeletal system morphogenesis (GO:0048705). This gene is essential during organogenesis and in particular during axial development in mice (Beck-Cormier et al., 2014). Several genes associated with other morphological traits in Maremmana likely went undetected by the present study either because of insufficient marker density or poor annotation of some of our candidate regions (e.g., the candidate region on the BTA12). The complete functional annotation of the bovine genome and the use of higher density SNP panels or sequencing data would be particularly relevant to improve our fine mapping results.

## CONCLUSION

In this study, signal consistency across three EHH-based methods was used to identify relevant genomic regions putatively under differential selection between Maremmana and a group of seven other Podolian-derived Italian breeds. Our results suggest that natural selection has shaped several immunity genes involved in both innate and adaptive response to the pathogens inhabiting Maremmana's natural habitat. Candidate genes associated with growth and meat traits appear to be targets of differential artificial selection between Maremmana and the other Podolian-derived Italian breeds. Information about the location of these regions can be used to detect variants that may underlie important traits of practical interest for the beef industry, and can help improving the competitiveness of the Maremmana cattle, through the implementation of genetic improvement programs for meat traits. Moreover, the candidate genes that confer a selective advantage for adaptative traits provide valuable information for future functional characterization as a starting point to identify causal genetic variants that control environmental adaptation traits for the genetic improvement of commonly used commercial breeds.

## DATA AVAILABILITY STATEMENT

The Illumina BeadChip 50 k SNP genotype data used and analyzed during the current study for the seven cattle breeds are deposited and available at https://www.animalgenome.org/repository/pub/UPIT2018.0607/. The genotyping data of the Maremmana breed presented in this article are not readily available because the raw datasets are property of the Consiglio per la Ricerca in Agricoltura e l'Analisi dell'Economia Agraria CREA and this information is commercially sensitive.

## AUTHOR CONTRIBUTIONS

SM conceived, planned, and coordinated the study. SB-J, GS, MB, and SM analyzed the data and performed the statistical analysis.
EC, RC, GC, and BP contributed with generation of data. FP and
BP acquired the funding. SB-J, EC, FP, and SM contributed to
the data interpretation. SB-J, GS, EC, RC, GC, MB, FP, BP, and SM
discussed the results and made suggestions and corrections. All authors
contributed to the article and approved the submitted version.

FUNDING

This research was financed by Livestock Environment Opendata
(LEO), PSRN 2014-2020 - Sottomisura: 16.2, Project number
PRJ-0185, CUP: J8418000090007.

REFERENCES

Ajmone-Marsan, P., Negrini, R., Milanesi, E., Bozzi, R., Nijman, I. J., Buntjer, J. B., et al. (2002). Genetic distances within and across cattle breeds as indicated by biallelic AFLP markers. *Anim. Genet.* 33, 280–286. doi: 10.1046/j.1365-2052.2002.00865.x

Bahbahani, H., Clifford, H., Wragg, D., Mbole-Kariuki, M. N., Van Tassell, C., Sonstegard, T., et al. (2015). Signatures of positive selection in East African Shorthorn Zebu: a genome-wide single nucleotide polymorphism analysis. *Sci. Rep.* 5, 1–13. doi: 10.1038/srep11729

Bal, M. S., Mahajan, V., Filia, G., Kaur, P., and Singh, A. (2016). Diagnosis and management of bovine babesiosis outbreaks in cattle in Punjab state. *Vet. World* 9, 1370–1374. doi: 10.14202/vetworld.2016.1370-1374

Band, G., Le, Q. S., Clarke, G. M., Kivinen, K., Hubbart, C., Jeffreys, A. E., et al. (2019). Insights into malaria susceptibility using genome-wide data on 17,000 individuals from Africa, Asia and Oceania. *Nat. Commun.* 10:5732. doi: 10.1038/s41467-019-13480-z

Beck-Cormier, S., Escande, M., Souilhol, C., Vandormael-Pournin, S., Sourice, S., Pilet, P., et al. (2014). Notchless is required for axial skeleton formation in mice. *PLoS One* 9:e98507. doi: 10.1371/journal.pone.0098507

Biscarini, F., Mastrangelo, S., Catillo, G., Senczuk, G., and Ciampolini, R. (2020). Insights into genetic diversity, runs of homozygosity and heterozygosity-rich regions in Maremmana semi-feral cattle using pedigree and genomic data. *Animals* 10:2285. doi: 10.3390/ani10122285

Boccolini, D., Toma, L., Di Luca, M., Severini, F., Cocchi, M., Bella, A., et al. (2012). Impact of environmental changes and human-related factors on the potential malaria vector, Anopheles branchiae (Diptera: Culicidae), in Maremma, Central Italy. *J. Med. Entomol.* 49, 833–842. doi: 10.1603/ME11252

Bomba, L., Nicolazzi, E. L., Milanesi, M., Negrini, R., Mancini, G., Biscarini, F., et al. (2015). Relative extended haplotype homozygosity signals across breeds reveal dairy and beef specific signatures of selection. *Genet. Sel. Evol.* 47:25. doi: 10.1186/s12711-015-0113-9

Bongiorni, S., Gruber, C. E. M., Chillemi, G., Bueno, S., Failla, S., Moioli, B., et al. (2016). Skeletal muscle transcriptional profiles in two Italian beef breeds, Chianina and Maremmana, reveal breed specific variation. *Mol. Biol. Rep.* 43, 253–268. doi: 10.1007/s11033-016-3957-3

Carmelo, V. A. O., and Kadarmideen, H. N. (2020). Genome regulation and gene interaction networks inferred from muscle transcriptome underlying feed efficiency in pigs. *Front. Genet.* 11:650. doi: 10.3389/fgene.2020.00650

Catillo, G., Moioli, B., Napolitano, F., and Steri, R. (2018). Identification of genomic regions harboring diversity between Holstein and two local endangered breeds, Modenese and Maremmana. *Livest. Sci.* 216, 75–83. doi: 10.1016/j.livsci.2018.07.011

Chen, W.C. (2011). Overlapping Codon Model, Phylogenetic Clustering, and Alternative Partial Expectation Conditional Maximization Algorithm. [Phd thesis]. Ames, Iowa State University.

Chen, Q., Huang, B., Zhan, J., Wang, J., Qu, K., Zhang, F., et al. (2020). Whole-genome analyses identify loci and selective signals associated with body size in cattle. *J. Anim. Sci.* 98:skaa068. doi: 10.1093/jas/skaa068

Ciecierska, A., Motyl, T., and Sadkowski, T. (2020). Transcriptomic profile of primary culture of skeletal muscle cells isolated from semitendinosus muscle of beef and dairy bulls. *Int. J. Mol. Sci.* 21:4794. doi: 10.3390/ijms21134794

ACKNOWLEDGMENTS

We thank BOVITA project for the agreement provided
among all the National Breeders Association and the
Universities involved.

SUPPLEMENTARY MATERIAL

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

Coble, D. J., Fleming, D., Persia, M. E., Ashwell, C. M., Rothschild, M. F., Schmidt, C. J., et al. (2014). RNA-seq analysis of broiler liver transcriptome reveals novel responses to high ambient temperature. *BMC Genomics* 15:1084. doi: 10.1186/1471-2164-15-1084

Consiglio Nazionale delle Ricerche (1983). Atlante Etnografico Delle Popolazioni Bovine Alleviate in Italia. Roma, Italy.

de Villemereuil, P., Frichot, É., Bazin, É., François, O., and Gaggiotti, O. E. (2014). Genome scan methods against more complex models: when and how much should we trust them? *Mol. Ecol.* 23, 2006–2019. doi: 10.1111/mec.12705

Di Lorenzo, P. D., Lancioni, H., Ceccobelli, S., Colli, L., Cardinali, I., Karsli, T., et al. (2018). Mitochondrial DNA variants of Podolian cattle breeds testify for a dual maternal origin. *PLoS One* 13:e0192567. doi: 10.1371/journal.pone.0192567

Donkervoort, S., Kutzner, C. E., Hu, Y., Lornage, X., Rendu, J., Stojkovic, T., et al. (2020). Pathogenic variants in the myosin chaperone UNC-45B cause progressive myopathy with eccentric cores. *Am. J. Hum. Genet.* 107, 1078–1095. doi: 10.1016/j.ajhg.2020.11.002

Du, X. H., Gan, Q. F., Yuan, Z. R., Gao, X., Zhang, L. P., Gao, H. J., et al. (2013). Polymorphism of MyoD1 and Myf6 genes and associations with carcass and meat quality traits in beef cattle. *Genet. Mol. Res.* 12, 6708–6717. doi: 10.4238/2013.December.13.4

Efergan, A., Azouz, N. P., Klein, O., Noguchi, K., Rothenberg, M. E., Fukuda, M., et al. (2016). Rab12 regulates retrograde transport of mast cell secretory granules by interacting with the RILP-dynein complex. *J. Immunol.* 196, 1091–1101. doi: 10.4049/jimmunol.1500731

Esche, C., Stellato, C., and Beck, L. A. (2005). Chemokines: key players in innate and adaptive immunity. *J. Investig. Dermatol.* 125, 615–628. doi: 10.1111/j.0022-202X.2005.23841.x

Fernandez, E. E., Oltjen, J. W., and Sainz, R. D. (2020). Mitochondrial abundance and function in muscle from beef steers with divergent residual feed intakes. *Animal* 14, 560–565. doi: 10.1017/S1751731119002209

Fijarczyk, A., and Babik, W. (2015). Detecting balancing selection in genomes: limits and prospects. *Mol. Ecol.* 24, 3529–3545. doi: 10.1111/mec.13226

Fioretti, M., Negrini, R., Biffani, S., Quaglia, A., Valentini, A., and Nardone, A. (2020). Demographic structure and population dynamics of Maremmana cattle local breed after 35 years of traditional selection. *Livest. Sci.* 232:103903. doi: 10.1016/j.livsci.2019.103903

Fleming, D. S., Weigend, S., Simianer, H., Weigend, A., Rothschild, M., Schmidt, C., et al. (2017). Genomic comparison of indigenous African and northern European chickens reveals putative mechanisms of stress tolerance related to environmental selection pressure. *G3: Genes, Genomes, Genetics* 7, 1525–1537. doi: 10.1534/g3.117.041228

Food and Agriculture Organization (1991). Animal Genetic Resources Information. Available at: http://www.fao.org/3/a-u6405t.pdf (Accessed January 14, 2021).

Food and Agriculture Organization (2002). Organic Agriculture and Genetic Resources for Food and Agriculture. Corporate Document Repository. Available at: http://www.fao.org/docrep/005/AC784E/ac784e-04.htm (Accessed January 14, 2021).

Gao, H., Antony, R., Srinivasan, R., Wu, P., Wang, X., and Li, Y. (2020). UCHL1 regulates oxidative activity in skeletal muscle. *PLoS One* 15:e0241716. doi: 10.1371/journal.pone.0241716

Gao, Y., Gautier, M., Ding, X., Zhang, H., Wang, Y., Wang, X., et al. (2017). Species composition and environmental adaptation of indigenous Chinese cattle. *Sci. Rep.* 7:16196. doi: 10.1038/s41598-017-16438-7
Gautier, M., and Vitalis, R. (2012). Rehh: an R package to detect footprints of selection in genome-wide SNP data from a haplotype structure. *Bioinformatics* 28, 1176-1177. doi: 10.1093/bioinformatics/bts115

Geserick, P., Kaiser, F., Klemm, U., Kaufmann, S. H. E., and Zerrahn, J. (2004). Modulation of T cell development and activation by novel members of the Schlafen (slfn) gene family harbouring an RNA helicase-like motif. *Int. Immunol.* 16, 1535-1548. doi: 10.1093/intimm/dxh155

Gharib, S. A., Johnston, L. K., Huizar, I., Birkland, T. P., Hanson, J., Wang, Y., et al. (2014). MMP28 promotes macrophage polarization toward M2 cells and augments pulmonary fibrosis. *J. Leukoc. Biol.* 95, 9-18. doi: 10.1189/jlb.1112587

Glass, E. J., Preston, P. M., Springbett, A., Craigmile, S., Kirvar, E., Wilkie, G., et al. (2005). Bos taurus and *Bos indicus* (Sahiwal) calves respond differently to infection with *Theileria annulata* and produce markedly different levels of acute phase proteins. *Int. J. Parasitol.* 35, 337-347. doi: 10.1016/j.ijpara.2004.12.006

Hay, E. H., and Roberts, A. (2018). Genome-wide association study for carcass traits in a composite beef cattle breed. *Livest. Sci.* 213, 35-43. doi: 10.1016/j.livsci.2018.04.018

Hu, Y., Su, J., Cheng, L., Lan, D., and Li, D. (2020). Pectoral muscle transcriptome analyses reveal high-altitude adaptations in Tibetan chickens. *Anim. Biol.* 70, 1-16. doi: 10.1163/15707563-bja10026

Jernigan, N. L., Paffett, M. L., Walker, B. R., and Resta, T. C. (2009). ASIC1 contributes to pulmonary vascular smooth muscle store-operated Ca²⁺ entry. *Am. J. Physiol. Lung Cell. Mol. Physiol.* 297, L271-L285. doi: 10.1152/ajplung.00020.2009

Kim, E.-S., Ros-Freixedes, R., Pena, R. N., Baas, T. J., Estany, J., and Rothschild, M. F. (2015). Identification of signatures of selection for intramuscular fat and backfat thickness in two Duroc populations. *J. Anim. Sci.* 93, 3292-3302. doi: 10.2527/jas.2015-8879

Kumar, A., Ashraf, S., Goud, T. S., Grewal, A., Singh, S. V., Yadav, B. R., et al. (2015). Expression profiling of major heat shock protein genes during different seasons in cattle (*Bos indicus*) and buffalo (*Bubalus bubalis*) under tropical climatic condition. *J. Therm. Biol.* 51, 55-64. doi: 10.1016/j.jtherbio.2015.03.006

Lei, H. (2019). Impact of Genetics on Meat Quality of Pigs and Beef Cattle. [Phd thesis]. Edmonton, Alberta: University of Alberta.

Li, Z., Chen, J., Wang, Z., Pan, Y., Wang, Q., Xu, N., et al. (2016). Detection of selection signatures of population-specific genomic regions selected during domestication process in Jinhua pigs. *Anim. Genet.* 47, 672-681. doi: 10.1111/age.12475

Li, M., Kao, E., Gao, X., Sandig, H., Limmer, K., Pavon-Eternod, M., et al. (2012). Codon-usage-based inhibition of HIV protein synthesis by human schlafen 11. *Nature* 491, 125-128. doi: 10.1038/nature11433

Li, M., Tian, S., Jin, L., Zhou, G., Li, Y., Zhang, Y., et al. (2013). Genomic analyses identify distinct patterns of selection in domesticated pigs and Tibetan wild boars. *Nat. Genet.* 45, 1431-1438. doi: 10.1038/ng.2811

Lin, Y.-H., Zhen, Y.-Y., Chien, K.-Y., Lee, I.-C., Lin, W.-C., Chen, M.-Y., et al. (2017). LIMCH1 regulates nonmuscle myosin-II activity and suppresses cell migration. *Mol. Biol. Cell* 28, 1054-1065. doi: 10.1091/mbc.E15-04218

Lonergan, E. H., Zhang, W., and Lonergan, S. M. (2010). Biochemistry of postmortem muscle — lessons on mechanisms of meat tenderization. *Meat Sci.* 86, 184-195. doi: 10.1016/j.meatsci.2010.05.004

Ma, Y.-F., Han, X.-M., Huang, C.-P., Zhong, L., Adeola, A. C., Irwin, D. M., et al. (2019). Population genomics analysis revealed origin and high-altitude adaptation of Tibetan pigs. *Sci. Rep.* 9:11463. doi: 10.1038/s41598-019-47711-6

Maretto, F., Ramljak, J., Sbarra, F., Penasa, M., Mantovani, R., Ivanković, A., et al. (2012). Genetic relationships among Italian and Croatian Podolian cattle breeds assessed by microsatellite markers. *Livest. Sci.* 150, 256-264. doi: 10.1016/j.livsci.2012.09.011

Mastrangelo, S., Ciani, E., Ajmone Marsan, P., Bagnato, A., Battaglini, L., Bozzi, R., et al. (2018). Conservation status and historical relatedness of Italian cattle breeds. *Genet. Sel. Evol.* 50:35. doi: 10.1186/s12711-018-0406-x

Mastrangelo, S., Jemaa, S. B., Ciani, E., Sottile, G., Moscarelli, A., Boussaha, M., et al. (2020a). Genome-wide detection of signatures of selection in three Valdostana cattle populations. *J. Anim. Breed. Genet.* 137, 609-621. doi: 10.1111/jbg.12476

Mastrangelo, S., Tolone, M., Ben Jemaa, S., Sottile, G., Di Gerlando, R., Cortés, O., et al. (2020b). Refining the genetic structure and relationships of European cattle breeds through meta-analysis of worldwide genomic SNP data, focusing on Italian cattle. *Sci. Rep.* 10:14522. doi: 10.1038/s41598-020-71375-2

Mohanarao, G., Mukherjee, A., Banerjee, D., Gohain, M., Dass, G., Brahma, B., et al. (2014). HSP70 family genes and HSP27 expression in response to heat and cold stress in vitro in peripheral blood mononuclear cells of goat (*Capra hircus*). *Small Rumin. Res.* 116, 94-99. doi: 10.1016/j.smallrumres.2013.10.014

Moioli, B., Napolitano, F., and Catillo, G. (2004). Genetic diversity between Piedmontese, Maremma, and Podolica cattle breeds. *J. Hered.* 95, 250-256. doi: 10.1093/jheresh/esh032

Newman, J. H., Holt, T. N., Hedges, L. K., Womack, B., Memon, S. S., Willers, E. D., et al. (2011). High-altitude pulmonary hypertension in cattle (brisket disease): candidate genes and gene expression profiling of peripheral blood mononuclear cells. *Pulm. Circ.* 1, 462-469. doi: 10.4103/2045-8932.93545

Ning, T., Xiao, H., Li, J., Hua, S., and Zhang, Y. P. (2010). Adaptive evolution of the mitochondrial ND6 gene in the domestic horse. *Genet. Mol. Res.* 9, 144-150. doi: 10.4238/vol9-1gmr705

Ojano-Dirain, C., Pumford, N. R., Iqbal, M., Wing, T., Cooper, M., and Bottje, W. G. (2005). Biochemical evaluation of mitochondrial respiratory chain in duodenum of low and high feed efficient broilers. *Poult. Sci.* 84, 1926-1934. doi: 10.1093/ps/84.12.1926

Online Mendelian Inheritance in Man (2003). OMIM (TM). Johns Hopkins University, Baltimore, MD. MIM Number: 601025. Available at: https://www.omim.org/entry/601025 (Accessed January 25, 2021).

Ouyang, H., Zhang, H., Li, W., Liang, S., Jebessa, E., Abdalla, B. A., et al. (2016). Identification, expression and variation of the GNPD2 gene, and its association with body weight and fatness traits in chicken. *PeerJ* 4:e2129. doi: 10.7717/peerj.2129

Paradis, E., Claude, J., and Strimmer, K. (2004). APE: analyses of phylogenetics and evolution in R language. *Bioinformatics* 20, 289-290. doi: 10.1093/bioinformatics/btg412

Pariset, L., Mariotti, M., Nardone, A., Soysal, M. I., Ozkan, E., Williams, J. L., et al. (2010). Relationships between Podolic cattle breeds assessed by single nucleotide polymorphisms (SNPs) genotyping. *J. Anim. Breed. Genet.* 127, 481-488. doi: 10.1111/j.1439-0388.2010.00868.x

Porter, V., Alderson, L., Hall, S. J. G., and Sponenberg, D. P. (2016). *Mason's World Encyclopedia of Livestock Breeds and Breeding*. Vol. 2. Oxfordshire: CABI.

Purcell, S., Neale, B., Todd-Brown, K., Thomas, L., Ferreira, M. A. R., Bender, D., et al. (2007). PLINK: A tool set for whole-genome association and population-based linkage analyses. *Am. J. Hum. Genet.* 81, 559-575. doi: 10.1086/519795

Qanbari, S., Gianola, D., Hayes, B., Schenkel, F., Miller, S., Moore, S., et al. (2011). Application of site and haplotype-frequency based approaches for detecting selection signatures in cattle. *BMC Genomics* 12:318. doi: 10.1186/1471-2164-12-318

Ramos, P. M., Li, C., Elzo, M. A., Wohlgemuth, S. E., and Scheffler, T. L. (2020). Mitochondrial oxygen consumption in early postmortem permeabilized skeletal muscle fibers is influenced by cattle breed. *J. Anim. Sci.* 98:skaa044. doi: 10.1093/jas/skaa044

Razavi, S. M., Nazifi, S., Hasani, S., and Rakhshandehroo, E. (2015). Bovine serum theileriosis: effects on the cardiovascular system on the basis of serum analysis. *Comp. Clin. Pathol.* 24, 29-33. doi: 10.1007/s00580-013-1850-1

Renström, F., Payne, F., Nordström, A., Brito, E. C., Rolandsson, O., Hallmans, G., et al. (2009). Replication and extension of genome-wide association study results for obesity in 4923 adults from northern Sweden. *Hum. Mol. Genet.* 18, 1489-1496. doi: 10.1093/hmg/ddp041

Ritchie, H. (2009). Breeds of Beef and Multi-Purpose Cattle: Acknowledgements. Available at: http://msue.anr.msu.edu/uploads/236/58553/breedsofbeefcattle_ritchie.pdf (Accessed January 25, 2021).

Romi, R., Boccolini, D., Vallorani, R., Severini, F., Toma, L., Cocchi, M., et al. (2012). Assessment of the risk of malaria re-introduction in the Maremma plain (Central Italy) using a multi-factorial approach. *Malar. J.* 11:98. doi: 10.1186/1475-2875-11-98

Rothhammer, S., Seichter, D., Förster, M., and Medugorac, I. (2013). A genome-wide scan for signatures of differential artificial selection in ten cattle breeds. *BMC Genomics* 14:908. doi: 10.1186/1471-2164-14-908

Ryu, J., and Lee, C. (2016). Genetic association of marbling score with intragenic nucleotide variants at selection signals of the bovine genome. *Animal* 10, 566-570. doi: 10.1017/S1751731115002633

Sabeti, P. C., Varilly, P., Fry, B., Lohmueller, J., Hostetter, E., Cotsapas, C., et al. (2007). Genome-wide detection and characterization of positive selection in human populations. *Nature* 449, 913-918. doi: 10.1038/nature06250
Scheet, P., and Stephens, M. (2006). A fast and flexible statistical model for large-scale population genotype data: applications to inferring missing genotypes and haplotypic phase. *Am. J. Hum. Genet.* 78, 629–644. doi: 10.1086/502802

Senczuk, G., Mastrangelo, S., Ajmone-Marsan, P., Becskei, Z., Colangelo, P., Colli, L., et al. (2020). Into the origin and diversification of Podolian cattle breeds by testing Neolithic scenarios of European colonization using genome-wide SNP data. *Authorea* [Preprint]. doi: 10.22541/au.160303335.52950496/v1

Shin, M. K., Shin, S. W., Jung, M., Park, H., Park, H. E., and Yoo, H. S. (2015). Host gene expression for *Mycobacterium avium subsp. paratuberculosis* infection in human THP-1 macrophages. *Pathog. Dis.* 73:ftv031. doi: 10.1093/femspd/ftv031

Signer-Hasler, H., Burren, A., Neuditschko, M., Frischknecht, M., Garrick, D., Stricker, C., et al. (2017). Population structure and genomic inbreeding in nine Swiss dairy cattle populations. *Genet. Sel. Evol.* 49, 1–13. doi: 10.1186/s12711-017-0358-6

Sorbolini, S., Marras, G., Gaspa, G., Dimauro, C., Cellesi, M., Valentini, A., et al. (2015). Detection of selection signatures in Piemontese and Marchigiana cattle, two breeds with similar production aptitudes but different selection histories. *Genet. Sel. Evol.* 47:52. doi: 10.1186/s12711-015-0128-2

Stairiker, C. J., Meurs, M.van, Leon, L. G., Brouwers-Haspels, A. A., Rijsbergen, L., Mueller, Y. M., et al. (2020). Heatr9 is an infection responsive gene that affects cytokine production in alveolar epithelial cells. *PLoS One* 15:e0236195. doi: 10.1371/journal.pone.0236195

Stella, A., Ajmone-Marsan, P., Lazzari, B., and Boettcher, P. (2010). Identification of selection signatures in cattle breeds selected for dairy production. *Genetics* 185, 1451–1461. doi: 10.1534/genetics.110.116111

Strillacci, M. G., Marelli, S. P., and Martinez-Velazquez, G. (2020). Hybrid versus autochthonous Turkey populations: homozygous genomic regions occurrences due to artificial and natural selection. *Animals* 10:1318. doi: 10.3390/ani10081318

Stroud, M. J., Nazgiewicz, A., McKenzie, E. A., Wang, Y., Kammerer, R. A., and Ballestrem, C. (2014). GAS2-like proteins mediate communication between microtubules and actin through interactions with end-binding proteins. *J. Cell Sci.* 127, 2672–2682. doi: 10.1242/jcs.140558

Tang, K., Thornton, K. R., and Stoneking, M. (2007). A new approach for using genome scans to detect recent positive selection in the human genome. *PLoS Biol.* 5:e171. doi: 10.1371/journal.pbio.0050171

Toma, L., Di Luca, M., Mancini, F., Severini, F., Mariano, C., Nicolai, G., et al. (2017). Molecular characterization of Babesia and Theileria species in ticks collected in the outskirt of Monte Romano, Lazio Region, Central Italy. *Ann. Ist. Super. Sanita* 53, 30–34. doi: 10.4415/ANN_17_01_07

Toma, L., Khoury, C., Bianchi, R., Severini, F., Mancini, F., Ciervo, A., et al. (2015). Preliminary investigation on tick fauna in the neighborhood of Tarquinia, Lazio, Italy. *Ann. Ist. Super. Sanita* 51, 67–70. doi: 10.4415/ANN_15_01_12

Velez-Irizarry, D., Casiro, S., Daza, K. R., Bates, R. O., Raney, N. E., Steibel, J. P., et al. (2019). Genetic control of longissimus dorsi muscle gene expression variation and joint analysis with phenotypic quantitative trait loci in pigs. *BMC Genomics* 20, 1–19. doi: 10.1186/s12864-018-5386-2

Vitti, J. J., Grossman, S. R., and Sabeti, P. C. (2013). Detecting natural selection in genomic data. *Annu. Rev. Genet.* 47, 97–120. doi: 10.1146/annurev-genet-111212-133526

Voight, B. F., Kudaravalli, S., Wen, X., and Pritchard, J. K. (2006). A map of recent positive selection in the human genome. *PLoS Biol.* 4:e72. doi: 10.1371/journal.pbio.0040072

Zhao, F., McParland, S., Kearney, F., Du, L., and Berry, D. P. (2015). Detection of selection signatures in dairy and beef cattle using high-density genomic information. *Genet. Sel. Evol.* 47:49. doi: 10.1186/s12711-015-0127-3

Zinovieva, N. A., Dotsev, A. V., Sermyagin, A. A., Deniskova, T. E., Abelmanova, A. S., Kharzinova, V. R., et al. (2020). Selection signatures in two oldest Russian native cattle breeds revealed using high-density single nucleotide polymorphism analysis. *PLoS One* 15:e0242200. doi: 10.1371/journal.pone.0242200

Zsolnai, A., Maróti-Agóts, Á., Kovács, A., Bálteanu, A. V., Kaltenecker, E., and Anton, I. (2020). Genetic position of Hungarian Grey among European cattle and identification of breed-specific markers. *Animal* 14, 1786–1792. doi: 10.1017/S1751731120000634

**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.

The handling editor declared a past collaboration with several of the authors EC, RC, FP, and SM.

Copyright © 2021 Ben-Jemaa, Senczuk, Ciani, Ciampolini, Catillo, Boussaha, Pilla, Portolano and Mastrangelo. 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.