Introduction

First documented in the 13^th Century on the Iberian Peninsula, the Lidia cattle breed is distinguished by selection for behavioral characteristics that enhance aggressiveness and for its use in civil and religious events¹. Various social and cultural phenomena involving bulls are currently included in what is known as “bullfighting”^2,3. Several countries consider different bullfighting traditions as practices that reinforce regional cultural identity^2,3; indeed, in Spain and Peru it has been designated an intangible cultural heritage⁴. The Lidia breed is characterized for having low genetic and ecological interchangeability^5,6.

The first documented bullfighting celebration in Mexico was held in 1523 using aggressive cattle brought mainly from the Navarra region in Spain, which is where the Casta Navarra breed originates⁷. It was not until the turn of the 20^th Century, however, that specialized breeding of Lidia began in Mexico with importation of a small number of animals from Spain by four breeding families: Llaguno, González, Barbabosa and Madrazo^8,9. Only genetic lines originating with the Llaguno and González families are still extant today^8,10.

The Llaguno family has been located largely in north-central Mexico. Under Antonio Llaguno the reproduction system was closed involving crosses only between Lidia animals directly linked to the original imported animals; in Mexican livestock terminology these are known as “Pure” animals. Julián Llaguno, brother of Antonio, followed a different breeding strategy, crossing Creole dams with Lidia sires of known Spanish origin; these are termed “Impure”^8,9. The González family, located in south-central Mexico, crosses imported Lidia breed animals with local cattle selected for aggressiveness⁸. With the purpose of breeding bulls for bullfights, in 1996 and 1997 a group of Mexican breeders imported animals from Spanish lines such as Domecq, Murube, Santa Coloma and Saltillo, among others; this strategy ended when livestock imports were prohibited for animal health reasons¹⁰.

The current Lidia breed population in Mexico is approximately 110,000 animals raised on a total of around 135,000 ha¹⁰. This breed is raised under extensive conditions, which favors conservation of endemic flora and fauna. Its central role in many local social traditions supports Mexico’s livestock economy while reinforcing regional cultural identity^7,10,11.

The genetic variability of the Lidia breed population in Mexico versus the original Spanish population has been analyzed using autosomal microsatellite markers, with differentiation between Spanish Lidia lines and the Llaguno and González family lines¹². These results were confirmed using molecular data produced with DNA chips for bi-allelic single nucleotide polymorphism (SNP) molecular markers. Clear genetic differentiation has also been reported between the Antonio Llaguno and González family lines^13,14, although these analyses did not include samples from the Julián Llaguno line or the lines imported in the late 20^th Century (Domecq, Santa Coloma, etc.).

The present study objective was to use SNP molecular markers to analyze the genomic diversity, population structure, endogamy levels and genetic relationships in representative populations of the Lidia breed in Mexico.

Material and methods

A total of 306 blood samples were randomly collected from animals belonging to 32 ranches in Mexico affiliated with the Union of Lidia Bull Breeders (Unión de Criadores de Toros de Lidia). The samples were classified into five lines based on the historical origins of each: Antonio Llaguno, Julián Llaguno, González, Domecq and Santa Coloma (Table 1). The samples were collected in tubes containing Magic Buffer^® preservative (Biogen Diagnostica, Spain) and kept at 15 °C until DNA extraction. Genomic DNA was extracted using a standard phenol/chloroform protocol¹⁵, and the samples were later genotyped with the 50K medium density SNP bovine chip (http://www.illumina.com).

Using the PLINK ver. 1.07 software¹⁶, the information was refined by excluding SNPs located on sex chromosomes, those exhibiting a minor (<0.01) allele frequency (MAF), those with <20% missing genotypes and those diverging from Hardy-Weinberg equilibrium (P<0.001). A total of 41,455 SNPs remained for analysis.

Again using PLINK¹⁶, analyses were done of three genetic diversity parameters: observed heterozygosity (Ho), expected heterozygosity (He) and the endogamy coefficient (F _IS ), estimated as 1-Ho/He. The F _ST coefficients were calculated using the ARLEQUIN ver. 3.0 software¹⁷. For each individual, the proportion of genetic origins identifiable using the Bayesian grouping algorithm was calculated with the ADMIXTURE software^18,19. Graphs were generated with the POPHELPER ver. 1.0.10 software²⁰.

A molecular analysis of variance (AMOVA) was run using a linear model to evaluate genetic variation between and within lines¹⁷. The analysis was done in hierarchical mode with three levels (between lines, between ranches in the same line and within ranches). The same software was used to calculate mean distance of the lines in terms of F _ST .

Individual runs of homozygosity (ROHs) were identified per individual²¹. This was done using PLINK with 30 SNP windows, allowing for <100 kb between two consecutive homozygous SNPs, less than two missing genotypes, one heterozygous and a 500 kbp minimum length. The average value per ranch and per line was then calculated.

Results and discussion

Genetic structure and population differentiation

The cross-validation error (CV) used in ADMIXTURE calculates values that decrease as the number of hypothetical ancestral populations (K) increases. When the CV value begins to increase it indicates the most probable hypothetical population prediction. Using the present data the most accurate prediction was identified at K = 5^18,19.

The average proportions of individuals in the ranches coincided with their assignment to each of the five ancestral populations (Figure 1). Each of the five lines largely corresponded to one of the five defined ancestral populations. Discrimination between the Julián Llaguno and Antonio Llaguno lines was less evident between some ranches in these lines, while it was greater between others (e.g., El Sauz and Valparaíso in Julián Llaguno, and Garfias, Los Cués and La Antigua in Antonio Llaguno). This analysis does not explain these differentiations between the Llaguno lines. Perhaps they result from variation in original genetic material since the Julián Llaguno line includes crosses between Creole dams and Spanish Lidia sires. Nonetheless, it is clear that both line (Antonio Llaguno and Julián Llaguno) mostly share common genetic origins. In contrast, the Gonzáles, Domecq and Santa Coloma lines are clearly genetically distinct.

Each vertical line represents an individual animal’s total genome. The proportion of each color (genetic group, K) in the vertical lines is the proportion of each of the five ancestral populations in an individual’s genome (K).

Figure 1 Cross-validation error analysis of hypothetical ancestral populations (K) using ADMIXTURE 

Identification of the ROHs per line produced statistics on the average number of runs or segments and average ROH length in each of the five lines (Figure 2). The number and length of ROHs in the different lines exhibited similar patterns. Greater number of segments and length of the ROH are correlated with recent consanguinity events²³.

Figure 2 Average number of runs of homozygosity (ROHs) and their average size (Mb) in the five studied lines of the Lidia breed in Mexico 

In the data for average ROH number and length by ranch (Table 3), length ranged from 5 (Jaral de Peñas) to 7.7 Mb (Boquilla del Carmen and El Sauz). This is consistent with the average F _IS values in which the highest values (>0.20) corresponded to Boquilla del Carmen and El Sauz while the lowest was for Jaral de Peñas (Table 1).

Table 3 Average number of runs of homozygosity (ROH) per ranch, including number of segments (NSEG) and average length (Mb) 

Line	Ranch	NSEG	Average length (Mb)
Julián Llaguno	Pozo Hondo	108	5.3
	Valparaiso	134	6.2
	El Sauz	111	7.7
	Caparica	114	6.5
	Pomedio	117	6.1
Antonio Llaguno	San Mateo	141	6.5
	Reyes Huerta	124	5.9
	Fernando de la Mora	87	5.8
	Garfias	139	7.1
	Los Cués	131	7.3
	La Antigua	144	6.9
	Xajay	115	6.6
	Teófilo Gómez	135	6.0
	Celia Barbabosa	126	5.8
	Boquilla del Cármen	132	7.7
	Fermín Rivera	131	6.1
	Corlomé	79	5.3
	Arroyo Zarco	124	6.3
	Marrón	105	6.2
	La Punta	104	6.1
	Pomedio	121	6.2
González	Tenexac	135	7.4
	Gonzalo Yturbe	110	6.3
	De Haro	113	5.6
	C.Castañeda	134	6.9
	Zacatepec	83	5.7
	Rancho Seco	95	5.2
	Pomedio	109	6.2
Domecq	La Joya	113	6.3
	Sta. Maria de Xalpa	92	5.7
	Jaral de Peñas	64	5.0
	José Julián Llaguno	72	5.1
	Torréon de Cañas	73	5.2
	Pomedio	85	5.5
Santa Coloma	Los Encinos	120	6.4
	San José	105	5.8
	Pomedio	112	6.1

When animals are isolated in relatively small populations the probability is greater that they inherit identical DNA segments that account for ROH²³. In previous studies, the high number and long length of ROHs have been associated with endogamy²³. This coincides with the results observed here for the five studied Lidia breed lines in Mexico, which had values greater than those reported in previous analyses of ROHs in native Spanish and American Creole breeds¹³. Both the ROH and F _IS values in the studied Mexican Lidia population reflect subdivision into lines and its consequences: reduction of effective sizes and higher consanguinity values²⁴. This subdivision is effective at preserving intrapopulation genetic variability²⁵,;however, as each subpopulation experiences genetic drift consanguinity will increase and genetic variability will decrease. Under these circumstances it is advisable to closely monitor degrees of endogamy.

Conclusions and implications

The genetic differentiation observed among the Mexican Lidia population, into lines and even between ranches, is due to the different genetic origins of some lines (i.e. Domecq and Santa Coloma) in conjunction with the genetic isolation maintained between the remaining lines (i.e. Antonio Llaguno, Julián Llaguno and González). Both the genetic structure and ROH analyses identified genetic isolation between the lines of the Mexican Lidia population, which has contributed to genomic variations when compared to European Lidia populations. The Mexican Lidia lines are clearly unique from their ancestral populations and quite differentiated amongst themselves.