Abstract
Domestic pigs were a key component of the Neolithic Revolution because of their great relevance to farming. Zoo-archaeological evidences suggest that Sus scrofa was domesticated in the Fertile Crescent about 10,500 years BP. From that moment, early Neolithic farmers spread domestic pigs westward into Europe. Yet, once domesticated, European pigs rapidly replaced pigs of Near Eastern origin throughout Europe. A temporal distribution change between European mitochondrial DNA haplotypes (A-side and C-side) also occurred: the A-side haplotype increased in domestic remains from the Neolithic to the Roman Age in Europe, at the expense of C-side individuals. This same pattern is absent in non-domestic settings. We jointly analyzed (modern) wild boar morphology and mitochondrial DNA, seeking out morphological differences between A- and C- side types. Our results show that A-side wild boars are significantly larger than C-sides, irrespective of sex, age, and reproductive stage. This suggests that the increased frequency of A-side individuals in domestic samples through time might be the direct result of active selection by early breeders for their fast growth rate.
Introduction
Animal domestication led to one of the most important socioeconomic transition in human history. During the Neolithic Revolution, the combination of herding and agriculture stimulated a rapid increase in human population size, which in turn spurred human migration from the centers of origin to their surroundings (Bocquet-Appel, 2011; Diamond, 2002; Diamond and Bellwood, 2003; Gignoux et al., 2011). This transition took place in a limited number of geographic regions across the globe, in the form of multiple, independent events (Purugganan and Fuller, 2009, 2011). Most livestock domestication occurred in Southwest Asia and East Asia around 11,000–8000 years BP and later in the Americas about 4500 years BP. The Fertile Crescent was the first center of domestic forms of wild sheep, goat, cow, and pig, as early as during the 10th millennium BC (Clutton-Brock, 1999; Zeder, 2006, 2012a, 2012b).
The wild boar (Sus scrofa) was an important prey to hunter-gatherers across wide areas of Eurasia until the early Holocene (Benecke, 1993; Pushkina and Raia, 2008), and currently, it is one of the most widespread wild mammals. Similarly, its domestic form, the pig, soon became of utmost economic importance. Domestic pigs are common in both archeological and paleontological contexts (Clutton-Brock, 1999; Masseti, 2002; Wilkens, 2003) and became economically relevant to humans from the late Pleistocene (Albarella et al., 2007; Masseti, 2007). Current pig varieties are the result of multiple domestication events that occurred in different world locations (Chen et al., 2007; Larson et al., 2005). Several studies employed a wide range of genetic markers to infer domestication and migration processes, and patterns of genetic diversification in both wild and domestic Sus scrofa (Aguilera-Reyes et al., 2006; Larson and Burger, 2013; Larson et al., 2014). Currently, mitochondrial DNA (mtDNA) analysis is the most common tool in domestication studies, for the rich data sets available and the advantages it confers to the analysis of ancient DNA (aDNA), including its fast evolutionary rate and the high number of copies of mtDNA in the eukaryote cell. Studies based on the mtDNA are thus optimally suited to analyze population divergence at the large scale, and for phylogeographic inference (Bruford et al., 2003; Larson and Burger, 2013; Larson et al., 2010; MacHugh and Bradley, 2001; MacHugh et al., 1997; Okumura et al., 2001).
The genetic diversity of Sus scrofa as deduced from mitochondrial control (D-loop) region analysis indicates the existence of three divergent lineages in Western Eurasia, geographically distributed in the Near East (NE), over most continental Europe (E1), and endemic to Italy (E2). The European clade (E1) further shows a geographical segregation into two haplogroups: A-side haplotypes are common in central Europe and Italy, while C-side haplotypes are mostly found in Iberia and Eastern Europe (Alex-andri et al., 2012; Alves et al., 2003, 2010; Giuffra et al., 2000; Larson et al., 2005; Scandura et al., 2011; Vilaça et al., 2014). The sharp geographical segregation of A-side versus C-side haplogroups in Europe changed over time, as a result of a post-Neolithic diffusion of pig lineages, and was probably absent before pig domestication (Larson et al., 2007; Scandura et al., 2011).
An extensive zoo-archeological record demonstrates that the wild boar was first domesticated in the Near East by at least 10,500 years BP (Conolly et al., 2011; Ervynck et al., 2001; Vigne et al., 2011). Recent molecular evidence suggests a second, independent domestication in Europe. By sequencing a short D-loop region from ancient remains, Larson et al. (2007) demonstrated that early Neolithic farmers introduced domestic pigs into Europe from the Near East by 6th millennium BC, through the so-called Danubian route, before spreading into the Paris basin as early as during 4th millennium BC. By ca. 3900 BC, European wild boars were locally domesticated and then spread throughout Europe and into the Near East.
By analyzing data in Larson et al. (2007), we noticed that the frequency of the A-side haplotype greatly increased in domestic contexts since the Neolithic, while the C-side haplotype concomitantly decreased in the same settings. Thus, A-side haplotype individuals seem to have been preferentially domesticated or more abundant either in domestic context (Appendix 1, available online). We postulated that this preference could reside in phenotypic differences between the two haplogroups. During domestication, wild individuals were selected to the advantage of humans (i.e. artificial selection). Deliberate human selection for desired traits (e.g. body size, meat or milk yield, coat quality) can also be recognized in the archaeological record (Zeder, 2006, Zeder et al., 2006). For example, chickens were selected to be larger, wild cattle (aurochs) to be smaller, and sheep to lose their bristly outer hair (the kemp) and not to shed their soft inner hair (the wool; Diamond, 2002). Grigson (1999) showed pig-raising for meat consumption at Windmill Hill prehistoric site. We hypothesized that A-side wild boar should either grow faster or be larger in order to be amenable to selection by early farmers. This would increase meat production rate to the advantage of breeders.
Since morphometric data on genotyped wild boar remains are exceedingly scarce, in order to verify this hypothesis, we analyzed the morphology of extant S. scrofa and its mtDNA control region. If the A-side wild boar is phenotypically different from C-side individuals, active selection by early breeders would then be proven, and hence explain the great preponderance of the A-side haplotype in domestic breeds.
Material and methods
Morphological measurements
To confront A- and C-side phenotypes, we analyzed 55 specimens taken from a culling program implemented to reduce agricultural damage caused by wild boar. The samples belong to 10 locations in the territory of ‘Cilento, Vallo di Diano, and Alburni National Park’ (PNCVD-181,000 ha). All the activities were developed from 2009 to 2012, in accordance to the Italian national laws on culling (157/92 and 394/91) and the approval of the executive board of PNCVD (10/2008). The sample includes 34 females and 21 males. For each individual, we analyzed a total of 15 morphological and reproductive traits: sex, coat color, age, body weight, number of fetuses (bearing in gravid females), total number of nipples, number of active nipples, body length (BL), chest circumference (CC), snout-ear tips length (SEL), limb length (LL), height at the shoulder (HS), tail length (TaL), height at the withers (HW), and ear length (EL; Mattioli and De Marinis, 2009; Šprem et al., 2011; see Appendix 2, available online for the linear measurements and Appendix 3, available online).
Age was assessed from tooth eruption stages and replacement patterns (Baubet et al., 1994; Matschke, 1967). Wild boars were classified into 2 distinct age classes according to Pedone et al. (1995): sub-adult (between 1 and 2 years old), and adult (older than 2 years). Sub-adults were then removed from the analyses. The total weight of the individual was measured by a digital dynamometer.
All linear body measurements were taken with a tape rule, according to the procedures for biometric detection of the Istituto Superiore per la Protezione e Ricerca Ambientale (ISPRA Institute; Mattioli and De Marinis, 2009). With reference to the reproductive state, all females were included into either one of these two categories: reproductive and non-reproductive. Referring to the first group, two additional types were defined: pregnant females bearing fetuses, and lactating females, depending on the count of active nipples and the presence of fetuses at dissection. In the case of pregnant females, the body weight considered was that of the mother minus the total weight of the unborn fetuses (FernáNdez-Llario et al., 1999, 2005). Out of a total of 34 females analyzed, 20 were pregnant.
DNA extraction, PCR, and sequencing
Muscle tissues collected upon culling were stored in 70% ethanol at −20°C. Genomic DNA was extracted using the DNeasy Blood and Tissue Kit (QIAGEN, USA) according to the manufacturer’s instructions. A fragment (638 base pair) of the mtDNA control (D-loop) region was amplified using the primers H16108 and L15433 (Watanobe et al., 2001). This region corresponds to positions 15,451–16,088 of the reference mtDNA genome (Ursing and Arnason, 1998, accession number: AJ002189). This region was also selected to maximize the size of possible alignments including already published GenBank sequences. Polymerase chain reactions (PCRs) were set up in 25 µL volumes containing: 25 ng of DNA template, 1 U of Taq DNA polymerase (QIAGEN, USA), 2.5 µL of 10X PCR buffer, 2.5 mM MgCl2, 0.25 mM of each deoxynucleotide triphosphate (dNTP), and 12.5 pmol of each primer. The amplification was performed in a Mastercycler (Eppendorf, Hamburg, Germany) under the following conditions: an initial DNA denaturation at 95°C for 3 min; 40 cycles at 94°C for 30 s, 62°C for 40 s, 72°C for 1 min, and a final elongation step at 72°C for 10 min. Amplifications of blank extractions and PCR blanks were performed in all experiments to monitor the contamination. PCR products of 5 µL were then separated by electrophoresis in a 1.5% agarose gel. PCR products were treated with Exonuclease I and Fast Alkaline Phosphatase (ThermoScientific, USA) in a Mastercycler at 37°C for 30 min and at 80°C for 15 min. The sequencing reaction was conducted in ABI 3100 automated DNA sequencer (Perkin-Elmer) using ABI Prism Big Dye Terminator v3.1 Cycle Sequencing kit (Applied Biosystems). The primers used for amplification were also used for sequencing, and each fragment was sequenced in both directions. Electropherograms were assembled to examine any base pair ambiguities and the new sequences aligned using software, Geneious version 5.4.3 (Biomatters, available from http://www.geneious.com/). DNA sequences were deposited in GenBank (http://www.ncbi.nlm.nih.gov/; see Appendix 3, available online).
Phylogenetic analysis of D-loop sequences
To ascribe culled individuals to any of the mtDNA haplogroups present in Europe, we first downloaded from GenBank 55 mtDNA control region sequences (see Appendix 4, available online) belonging to European Sus scrofa (Alves et al., 2003; Giuffra et al., 2000; Gongora et al., 2003; Kim et al., 2002; Larson et al., 2005; Randi et al., 2002; Watanobe et al., 1999, 2001). Bayesian phylogenetic analyses were carried out in MrBayes 3.2 (Ronquist and Huelsenbeck, 2003). The best fit model of nucleotide substitution was estimated using a hierarchical likelihood ratio test (hLRT) as implemented in jMODELTEST 3.7 (Guindon and Gascuel, 2003; Posada and Crandall, 1998). The best model resulted the HKY model with inv-gamma distributed (I+G) rate variation across the sites. The best fit parameters were included as priors for MrBayes 3.1.2 and three Markov Chain Monte Carlo (MCMC) searches were run for 15 million generations. For each MCMC, a tree was sampled every 100 generations. The first 10% of the trees were discarded as burn-in. A consensus tree from the retained trees was computed with MrBayes. Sus celebensis D-loop sequence was included as outgroup (GenBank accession number: AY884738). Sequences belonging to European clade were selected to create a new sequence data set and collapsed in haplotypes using the DnaSP version 5 software (Librado and Rozas, 2009). In order to identify the A-side and C-side haplotypes, a Median Joining network of haplotypes was constructed using only the European samples with the software Network 4.5.1.6 (Fluxus Technology, http://www.fluxusengineering.com/).
Analysis of variance
We analyzed adult specimens metric measurements to investigate upon phenotypic differences between the two E1 haplotypes. Some variables (body weight, body length, chest circumference, and limb lengths) represent the individual ‘size’ while others (tail length, ear length, and snout-ear tips length) are more representative of the body ‘shape’. We conducted a multivariate analysis of variance (MANOVA) including all of the variables. Then, we tested individual variables by means of univariate analyses of variance (ANOVAs). In each analysis of variance, test design, sex, haplotype appurtenance, and the interaction between the two were taken as factors. We considered the distribution of haplotypes in the 20 pregnant individuals, to exclude the influence of reproductive status on morphological differences between the haplotypes.
Principal component analysis
To reduce data dimensionality and account for normality, we further performed a Principal Component Analysis (PCA) on morphological data (Meloro et al., 2013). PCA scores were then used to perform MANOVA taking haplogroup assignment, sex, and their interaction as fixed factors.
Results
Identification of mtDNA haplotypes
We sequenced the mtDNA control region (638bp) in 55 living wild boars to identify the mitochondrial clade they belong to. The new sequences were aligned with 55 additional sequences representative of the current genetic diversity of Sus scrofa in Western Eurasia. Of the 55culled wild boar sequences, 52 sequences belong to the European E1 clade and 3 samples belong to the endemic Italian E2 clade (see Appendix 5, available online). E2 individuals were excluded from the analysis of haplotype identification. All the European sequences were selected for a new alignment and collapsed into distinct haplotypes (see Appendix 6, available online). A Median Joining network was created to discriminate European A-side from C-side haplotypes, using as a reference two haplotypes representing the two core lineages, separated by a transversion (Larson et al., 2005; see Appendix 7, available online). The samples were then assigned to either A- or C-side haplotypes. In 52 sequences of wild boar, 44% corresponded to the A-side haplotype (N = 23) and 56% to the C-side (N = 29).
Analysis of variance
MANOVA gave non-significant differences between A- and C-side individuals (factor haplotype: F = 1.546, p = 0.165; factor sex: F = 3.428, p = 0.003; sex/haplotype interaction: F = 2.084, p = 0.054). Yet, almost all of the ANOVA tests on ‘size’ variables indicated significant differences irrespective of sex. Most notably, A-side haplotype individuals were larger (mean body weight difference = 7.63 kg, p = 0.028), longer (mean body length difference = 10.12 cm, p = 0.027), and had a wider chest (mean chest difference = 6.340, p = 0.028). Interestingly, the two phenotypes did not differ in height at withers. Yet, we found positive allometry in both tail and ear lengths (Table 1). Since A-type wild boars are distinctively larger, their proportionally longer ears and tails (a phenotype that is recurrent among domestic breeds) would make them easy to recognize. Pregnancy status did not affect the results.
ANOVA test results for individual variables.
The PCA on logged data suggested minor differences between the two clades to occur (see Appendix 8, available online). In keeping with this, MANOVA conducted on PCA scores gave non-significant results (sex: Pillai trace = 0.169, F = 1.488, p = 0.204; haplogroup: Pillai trace = 0.162, F = 1.415, p = 0.230; interaction: Pillai trace = 0.154, F = 1.339, p = 0.261). Marginal differences are nonetheless linked to Principal Component 1 (which explains 41.2% of total variance), where scores differ per haplogroup (F = 4.615, p = 0.036). This component is mainly loaded on body size (r = −0.696) and height at the shoulder (r = −0.473). No additional significant differences in PCA scores are noted on other axes.
Discussion
The analysis of ancient DNA from pig remains revealed a complex temporal and geographical pattern of changes in pig genetic signature in Europe. Early Neolithic farmers spread domesticated pigs westward into Europe. Once domesticated, European mtDNA lineages replaced Near Eastern mtDNA signatures to spread throughout Europe. Regardless of the specific cause or progression of European pig domestication, European wild boar rapidly became the predominant lineage within the domestic forms. The ancient distribution of mitochondrial haplotypes revealed that during the wild boar domestication process (that is from the Neolithic onwards), the relative frequency of European A-side haplotype individuals increased at the expense of C-side individuals in domestic contexts (see data in Larson et al., 2007 and Appendix 1, available online). However, the reason of this substitution is unknown, which might be explained by morphological features good for human. Unfortunately, because of the nature of the fossil record, most of the ancient remains from hunting or domestic context are often fragmentary and it is not possible to deduce morphological traits such as weight, total length, age, or sexual dimorphism from these remains (Vigne et al., 2011, see also supplementary information in Larson et al., 2007). In order to analyze the wild reservoir on which artificial selection acted during the domestication, we were thus forced to study living boars, assuming that their morphology mirrors the past (a usual assumption in paleontological studies known as the actualism principle). We analyzed a range of biometric traits in current wild boar individuals belonging to either of the two haplotype groups. We showed that individuals belonging to the A-side haplotype are statistically larger, irrespective of sex, age, and reproductive status. Our hypothesis does not exclude the possibility that morphological differences between the two haplotypes partially depend on hybridization with domestic pigs nowadays (although we exclusively sampled wild individuals), even if the fixation of phenotypic traits of domestic individuals in wild populations happens to be exceedingly rare, and we see no clear reason why hybridization should occur in one haplotype more than on the other. Whatever the case, the presence of ‘domestic’ traits (such as faster growth rate) in wild individuals remains improbable since artificial selection produces phenotypes that are good to humans, not to survival in the wild (Zeder, 2012a, 2012b). Taken these caveats in mind, our results support the notion that early breeders did preferentially select larger individuals (hence preeminently A-side animals) to convert to domestication, and the signature is appreciable even today in wild populations. This preference is most likely linked to meat production. Meat consumption has a long history in human evolution, likely going back to the earliest known human-like ancestor living 5–7 million years ago (Lee-Thorp, 2002; Lee-Thorp et al., 1994; Sillen, 1992). It is well-known that most hunters focused on larger prey to maximize the energetic return of hunting (Zeder, 2006). Size was consciously increased in small domestic species, in order to the meet the demand for increased meat consumption (Clutton-Brock, 1999). It is therefore suggested here that the frequency change of the A-side at the expense of C-side in the domesticated pigs may have been the result of active selection to increase the economic performance of pig domestication. Early breeder had no hint about the genetics of the wild boars they attempted domesticating. Yet, we found A-type wild boars have distinctively longer ears (all else being equal), a phenotype that is recurrent among domestic breeds and easy to recognize.
Although we analyzed the relationship between morphology and phenotype in 3 samples only of wild boars belonging to the Italian clade, it is evident that these appear to be smaller than European clades representatives. The Italian clade was found only in the Peninsula and in Sardinia, probably because of a stronger isolation of the Italian populations (Vilaça et al., 2014). Further analysis is necessary to strengthen our contention that the Italian haplotypes are uncommon in domestic context because of their small size, which still persists today.
Despite the restricted area of the sampling we focused upon, the statistically significant evidence about the morphology of mtDNA haplotypes in our data can be interpreted as a valid tool to correlate a mitochondrial genetic marker and the body size, in the specific case of wild boar. The comparative analysis of mtDNA control region and size-related traits have never been used as a method to discriminate the morphological differences in wild species and could be applied in the future for other species. Because of fragmentary nature of ancient samples, this procedure can be applied with many advantages to studies which infer about past phenomena involving both the spread of domesticated species and selection of particular traits.
Footnotes
Acknowledgements
We are grateful to the ‘Cilento, Vallo di Diano, and Alburni National Park’ personnel for granting permits to conduct this study and for their kind assistance during the culling program.
Funding
The author(s) received no financial support for the research, authorship, and/or publication of this article.