The first Brazilian bovine breed: structure and genetic diversity of the Curraleiro Pé-duro

Sampling and collection of biological material

The sampling of the specimens used in this study was designed to avoid the Wahlund effect (Meeûs, 2018), and consisted of 474 individuals, representing approximately 10% of the known population of Curraleiro Pé-duro cattle.

Hair follicle samples from males and females were taken from various categories (calves, yearlings, cows and bulls), with the adults registered as base herd with the Associação Brasileiro de Criadores de Curraleiro Pé-duro.

The database is available at https://doi.org/10.5281/zenodo.7115391

The three herds visited (herein referred to as farms or subpopulations “A”, “B”, and “C”) are participants in the genetic improvement program of the breed, conducted by ABCPD, the Federal University of Piauí (FUP), and the State Univerity of Piauí (UESPI). This research project has been approved by the Research Ethics Committee of the Federal University of Piauí under #683/21 CEUA-UFPI. Farm “A” is located in the Carnaubais region of Piauí (−4.652482218991185, −42.05797250656697), Farm “B” in the Cocais region of Maranhão (−5.231024540320055, −44.47868128067367) and Farm “C” in the Sambito Valley of Piauí (−6.082350629527833, −42.24258929956412). Each farm represents a subpopulation, and exchanges of a small number of individuals occur irregularly and sporadically among farms, during the process of commercialization of animals.

Genotyping with microsatellite markers

PCR products were obtained with primers marked with fluorophores and, after purification, were examined with a Thermo Fisher ABI 3730 DNA analyzer (Thermo Fisher, Waltham, MA, USA). Fragment size was determined with the GeneMapper^® program, generating a file with the alleles present in each animal, for each marker tested.

The genotyping of the animals was outsourced to Laboratório Raça (Goiânia, Brazil), accredited by the Ministry of Agriculture, Livestock and Supply, Brazil (MAPA).

Seventeen microsatellite markers were used (Table 1) following MAPA #45 of 12/15/2017 (MAPA BRASIL, 2017), and in agreement with the results of comparative tests carried out by the International Society of Animal Genetics (ISAG). These markers are suitable for the investigation of the genetic structure, and the optimization of the findings aimed at the applicability of individual results in the identification of kinship, useful in the structuring of the breed that currently has only a base herd.

Data analysis

Micro-Checker 2.2.3 software (Van Oosterhout et al., 2004) was used to check for errors in the genotyping and tabulation of results, including the verification of possible monomorphic alleles, alleles outside the expected size range and the presence of stutter bands. Following this quality control, the data were analyzed based on the calculation of PIC, Hardy-Weinberg equilibrium test with Bonferroni correction at 5% significance level; and the estimation of null allele frequency performed with Cervus 3.0.7 software (Kalinowski, Taper & Marshall, 2007).

The number of alleles (Na), effective alleles (Ne), and expected (He) and observed (Ho) heterozygosity were estimated with the GenAIEx 6.5 software (Peakall & Smouse, 2012).

Wright’s F statistics (Wright, 1949) were carried out using the Gene-Pop Package version 4.7.5 (Rousset, 2022) of R 4.2 (R Core Team, 2022).

To examine differences among individuals and populations, the analysis of molecular variance (AMOVA) by Excoffier, Smouse & Quattro (1992) was performed based on the allele frequency of haplotypes with a nonparametric statistical test of permutations (Excoffier, Smouse & Quattro, 1992), using the Arlequin 3.5.2.2 software (Excoffier & Lischer, 2010) configured to carry out 10 thousand permutations using ΦST values.

Genetic differences among populations were examined as a function of the geographic distances among them using the Mantel test. Geographic distances among populations (Table 2) were log transformed to linearize the matrix of genetic distance and distance among populations (Smith & Weissman, 2020). In this test, the Slatikin’s genetic distance used was based on 10 thousand permutations in Arlequin 3.5.2.2 software (Excoffier & Lischer, 2010).

Structure 2.3.4 (Pritchard, Stephens & Donnelly, 2000) was used to perform a cluster analysis based on Bayesian statistics using the Monte Carlo Markov Chains method (MCMC) according to the similarity of genotypes, evaluating number of possible genetic groups (K) ranging from 1 to 10. The program was set to run 500,000 MCMC samples with 100,000 burn-in iterations (400,000 for inference) and for each run as indicated by Gilbert et al. (2012) on 20 repetitions.

Because animals are exchanged among properties (subpopulations), the program was adjusted to run the admixture model (Admixture) that aims to infer a single Alpha and use it for all populations. However, it was tested with an alpha estimate for each population (farm) (Data S1).

The results were analyzed using the main interpretive methodology by Evanno, Regnaut & Goudet (2005), and the most recently recommended method by Puechmaille (2016).

The graphical representation of the results was carried out using the Structure Selector software (Li & Liu, 2018).

Analysis of STR usability and data integrity

Seventeen microsatellite loci showed high level of polymorphism to evaluate the genetic diversity of the sampled populations (Table 3), with a mean PIC of 0.6184, considered as very informative (PIC >0.50) by the classification of Botstein et al. (1980).

Indications of null alleles were identified in only a portion of the loci used and in low incidence ($\bar{y}=0.016$). Therefore, most were real alleles, with a low incidence of null alleles (An). The highest incidence was in TGLA126, with an estimate of zero or very close to zero for others. This is an indication of the reliability of the laboratory analysis, ruling out the possibility that the genetic differentiation measures were overestimated. Evidence of null alleles was found at only a portion of the loci used, at low incidence.

Genetic variation of subpopulations

The estimated expected heterozygosity (He) was 0.74 (Table 4).

Based on the number of heterozygous individuals, the observed heterozygosity (Ho) was estimated, and was observed to be same as He.

Genetic differentiation among animals of the Curraleiro Pé-duro breed

The fixation indices proposed by Wright (Wright, 1949) are shown in Table 5.

Because F_IS is the product of the difference between He and Ho in relation to He, it indicates whether there are changes in population structure. Using a population in Hardy-Weiberg equilibrium (E-HW) as a model, the expected heterozygosity of the subpopulations (Farms A, B, C) was higher than expected in a theoretical population (Table 5). The low values obtained suggest a population close to equilibrium.

Genetic structure of the sampled population

The molecular analysis of variance (AMOVA) revealed the degree of genetic structure of the subpopulations sampled (Table 6).

The 98.5% of the variation arising from individuals within populations suggests that there is greater variation within the herd than among herd. Thus herds are relatively similar, with greater variation occurring among the individuals that make up each herd. This condition mischaracterizes the degree of genetic structuring among sampled population.

In general, the population analysis revealed a low structure in each sampled location (Table 7), but the F_ST values found were significant (P-value <0.05). Despite being low, the value observed for the population of Maranhão was higher than that of the populations of Piauí compared between them (Table 7). However, genetic distance might not be associated with geographic distance because the results of the Mantel test were not significant (Table 8) (p-value > 0.05.

F_ST ranging from 0.00723 and 0.03198); p-values <0.05)

Regarding differences among populations, despite AMOVA showing that the difference among herds is much smaller compared to intra-herd differences, there are indications that the genetic distance between the herd of Maranhão and the herds of Piauí is greater, based on Nei’s distance, statistically significant by F_ST (Table 7).

Genetic distance is not associated with geographic distance according to the Mantel test (Table 8).

Although almost the entire genetic distance among groups is explained by geographic distance, it is not statistically significant (p > 0.05). It should be pointed out that there are no other variables that notoriously affect the distance among groups (subpopulations). Therefore, there is no evidence of isolation by geographic distance.

Genetic differences among individuals

For the clustering of individuals according to similarity of genotypes, simulations of groupings were carried out with the Structure software (Pritchard, Stephens & Donnelly, 2000), testing variations in configurations, including disregarding the possibility of permutation among populations (farms) given that the commercialization of cattle among properties varies. However, convergence was observed when using independent alphas (CLUMPAK in supplements).

No single cluster was able to perfectly describe all the variability found in the sampled subpopulations. The simulation data analyzed by the method described by Evanno, Regnaut & Goudet (2005), which is based on Delta K values, indicated that the best grouping of genotypes occurs in two clusters (Fig. 1).

Gilbert et al. (2012) proposed the use of other methods of interpretation of structure simulations, especially when obtaining K = 2 with microsatellite data, such as expanding the methodologies for analyzing the results of the simulations, in order to avoid spurious clusters that erroneously project the value of K.

The values shown in Table 9 are in agreement with the number of clusters obtained with the method proposed by Evanno, Regnaut & Goudet (2005).

Therefore, in both methods, K = 2 was the appropriate cluster number to describe the number of subpopulations. The 10 simulations converged in the 20 trials performed with K = 2 (see supplement) indicating the same cluster, following Kopelman et al. (2015) for clusters that showed greater convergence.

Considering two main groups, the participation of each individual (vertical line) within each rectangle was evaluated into subgroupings in Fig. 2, which makes up the individualized grouping for each ΔK.

The genetic compositions of each of the subpopulations within each cluster is shown in Table 10, detailing the information in the above graph.