Analysis of genetic distance between Peruvian Alpaca (Vicugna Pacos) showing two distinct fleece phenotypes, Suri and Huacaya, by means of microsatellite markers

Abstract

Two coat phenotypes exist in Alpaca, Huacaya and Suri. The two coats show different fleece structure, textile characteristics and prices on the market. Although present scientific knowledge suggests a simple genetic model of inheritance, there is a tendency to manage and consider the two phenotypes as two different breeds. A 13 microsatellite panel was used in this study to assess genetic distance between Suri and Huacaya alpacas in a sample of non-related animals from two phenotypically pure flocks at the Illpa-Puno experimental station in Quimsachata, Peru. The animals are part of a germplasm established approximately 20 years ago and have been bred separately according to their coat type since then. Genetic variability parameters were also calculated. The data were statistically analyzed using the software Genalex 6.3, Phylip 3.69 and Fstat 2.9.3.2. The sample was tested for Hardy-Weinberg equilibrium (HWE) and after strict Bonferroni correction only one locus (LCA37) showed deviation from equilibrium (P<0.05). Linkage disequilibrium (LD) was also tested and 9 loci associations showed significant disequilibrium. Observed heterozygosis (Ho= 0.766; SE=0.044), expected heterozygosis (He=0.769; SE=0.033), number of alleles (Na=9.667, SE=0.772) and Fixation index (F=0.004; SE=0.036) are comparable to data from previous studies. Measures of genetic distance were 0.06 for Nei’s and 0.03 for Cavalli-Sforza’s. The analysis of molecular variance reported no existing variance between populations.

Considering the origin of the animals, their post domestication evolution and the reproductive practices in place, the results do not show genetic differentiation between the two populations for the studied loci.

Key Words:

• Alpaca
• Vicugna Pacos
• Microsatellite
• Genetic distance
• Genetic diversity.

Introduction

Two coat phenotypes exist in Alpaca, Huacaya and Suri. They differ phenotypically in terms of fleece structure. The first is a classically built fleece: compact, bulky and with high fibre crimp, similar to a merino fleece. The second is a fleece with longer and lustrous fibre organized in defined, hanging locks, more similar to a Lincoln or an Angora type of fleece.

These two phenotypes also differ with regard to fibre structure and characteristics (Antonini et al., 2001; Renieri et al., 2004; Frank et al., 2006). In terms of demography, in the Peruvian population there is a predominance of the Huacaya type which represents 90% of the alpaca fleece processed in Peru (Hoffman and Fowler, 1995). This occurs despite the fact that the Suri trait seems to segregate in a way closely similar to a single dominant gene or a haplotype (Ponzoni et al., 1997).

A different distribution of phenotypes can be observed in North America, Europe and Australia, where growing interest of both the textile industry and breeder associations for the Suri fibre has led to increasing Suri/Huacaya ratios and often to a paid premium for Suri fibre in the same fineness range as Huacaya (McGregor, 2006). Over the years this trend has generated two different schools of thought about the way the two phenotypes should be bred, classified and genetically managed. This is mainly true at the level of farming and breeder associations in North America and in Australia who often recommend considering and managing the two phenotypes as two separate breeds in order to keep the Suri line pure (Holt, 2005; Walker, 2006; Baychelier, 2002). This is not supported by present scientific knowledge, which suggests that the phenotype fleece type is a qualitative trait determined by 1 or 2 loci (Ponzoni et al., 1997; Renieri et al., 2009b; Presciuttini et al., 2010).

Post domestication evolution and secondary breed structure

Wheeler et al. (2006) recently confirmed the unanimous agreement that alpacas were domesticated from vicuñas (Vicugna vicugna mensalis) but the species’ post domestication evolution has been strongly affected by historical events. In fact, its evolution is divided into three main periods: i) pre-Conquest (from the domestication to the arrival of the conquerors); ii) the great Conquest crisis (lasting only a few decades, but with devastating consequences); and iii) the subsequent development (new expansion of domesticated camelids to the present day). The potential existence of breeds in the first period is difficult to study and remains unknown; however, it is important to stress that the Conquest caused a drastic decline in population size (90% according to Flores Ochoa, 1982), the geographic marginalization of the animals and the loss of any reproductive barrier. The complete disruption of the breeding system led to the mixing of species because four species co-existed within the same area: llamas, alpacas, vicuna and guanaco. The genetic consequences of these events were drastic and permanently changed the animals in comparison to the preceding specialization of pre-Conquest times. In this context, natural selection regained a dominant role over artificial selection and the specialization that had been established was lost. Since the Conquest, there has been a steady and slow demographic increase. At present, the demographic structure of the Andean alpaca population seems to be based on primary populations, also called primitive breeds, in which natural selection is favored over artificial selection. The mating system is panmictic and the phenotypic variation inside the flocks is high (Mason 1973; Denis 1982). For the moment, secondary standardized breeds do not seem to exist, apart from some experiences on very large private farms (Pacomarca, Mallkini). Nevertheless, the Andean alpaca population is subject to substantial genetic erosion through the diffusion of the full white phenotype. Suri is considered a phenotype with reduced fitness; for this reason, it is generally reared separately from Huacaya (Renieri et al., 2009a).

Population genetics by means of microsatellite analysis

A 10 microsatellite markers panel has been used by the ARI (Alpaca Registry Inc.) since 1998, mainly for parentage verification, and has found other applications throughout the years in other fields such as genome mapping, population structure and comparable genome analysis, as reviewed by Munyard et al. (2009). The continuous and ongoing effort in sequencing the alpaca genome will rapidly lead to a much larger set of markers, as demonstrated by a number of studies published in the past decade (Obreque et al., 1999; Penedo et al., 1999; McPartlan et al., 1998). Reed and Chaves (2008) report an additional 1516 potential loci by blasting bos taurus SSRs and Munyard et al. (2009) have recently found a set of 9 tetranucleotide markers. Some of these markers have already been used to calculate genetic distances among different species of South American Camelids (Wheeler et al., 2006; Bustamante et al., 2002). Since Goldstein et al. (1995) evaluated their use for the calculation of genetic distances, they have been used in an increasingly large number of species.

Given the unusual post domestication evolutionary history of the species and the increasing interest in Suri fibre, the aim of this study was, therefore, to use a microsatellite panel to study the genetic distance between Suri and Huacaya alpacas and to assess the amplitude of genetic variability in the Peruvian alpaca population. The studied population belongs to a germplasm established approximately 20 years ago in Quimsachata (Puno province, Peru) and has been bred separately according to their coat type since then.

Materials and methods

Material collection and sample structure

Alpaca blood samples were collected in spring 2008 within a larger sampling program including skin biopsies and fibre. Animals were kept and managed in the Illpa-Puno Experimental Station in Quimsachata (Puno province, Peru) at an altitude of approximately 4200 m a.s.l. The animals belonged to two geographically separate and phenotypically pure flocks (one Suri and one Huacaya) that have been managed and bred separately since the alpaca germplasm was created at the experimental station 20 years ago. It is important to emphasize that common breeding practices in local Peruvian communities involve rearing Suri and Huacaya animals separately.

Blood samples were taken from a subset of 65 non-related animals selected for the microsatellite analysis. The sample was structured in order to include an equal number of Suri and Huacaya individuals and to respect the sex ratio as much as possible. It included all the available males (n=15) and a subset of females (n=50): in total 32 Huacaya (7 males and 25 females) and 33 Suri (8 males and 25 females).

For convenience and due to the lack of basic facilities, blood was collected by spotting a total of 100 µL of blood on Whatman FTA Nucleic Acid Collection cards (# WB120205). Possibly due to the non-sterile and difficult conditions of facilities on the Peruvian plateau, not all samples allowed a sufficient quality/quantity of DNA to be amplified for all microsatellites, and in order to minimize the number of missing data, the dataset was rearranged to obtain a final number of 49 individuals: 10 males (5 Huacaya and 5 Suri) and 39 females (19 Huacaya and 20 Suri), and 13 microsatellites; LCA 19 was not included in the analysis.

Microsatellite amplification details

All samples were processed in Italy and genomic amplification was carried out by LGS genetic laboratories (Cremona, Italy). The markers, dyes utilized, primer sequences, multiplex associations and allele sizes for the 14 microsatellites are shown in Table 1. The panel has been optimized to be amplified in two multiplex reactions of 7 primer pairs each.

Table 1 Markers and multiplex reaction data.

Marker	Fragment length	Alleles	Dye	5’-3’ Primer pair sequences
LCA 19	80-122	17	Vic	taagtccagccccacactcaggtgaaggggcttgatcttc
LCA 94	187-213	9	Pet	gtccattcatccagcacaggacatttggcaatctctggagaa
YWLL 44	84-136	18	Ned	ctcaacaatgctagaccttgggagaacacaggctggtgaata
YWLL 36	136-176	17	Vic	agtcttggtgtggtggtagaatgccaggatactgacagtgat
YWLL 43	128-164	10	Pet	Atacctctcttgctctctctccctctacaaccatgttagcca
YWLL 29	210-232	9	Fam	gaaggcaggagaaaaggtagcagaggcttaataacttgcag
LCA 37	124-174	19	Fam	Aaacctaattacctcccccaccatgtagttgcaggacacg
LCA 5	178-218	13	Vic	Gtggtttttgcccaagctcacctccagtctggggatttc
LCA 8	211-261	14	Pet	gctgaaccacaatgcaaagaaatgcagatgtgcctcagtt
LCA 65	159-193	14	Fam	Tttttcccctgtggttgaataactcagctgttgtcagggg
LCA 66	216-266	24	Ned	gtgcagcgtccaaatagtcaccagcatcgtccagtattca
YWLL 40	176-190	7	Ned	cacatgaccatgtccccttatccagtgacagtgtgactaaga
LCA 99	263-297	11	Vic	caggtatcaggagacgggctagcatttatcaaggaacaccagc
YWLL 46	87-115	5	Fam	aagcagagtgatttaaccgtgggatgactaagactgctctga

List of markers, allele fragment lengths, number of alleles, dyes used and respective primer pairs in the multiplex reactions.

Amplifications were carried out in 20 µL reactions with final concentrations of 0.2 µM of each primer, 1X Ampli-Taq Gold Buffer, 1.5 mM MgCl₂, 200 µM of each dNTP and a total of 2 units of Ampli-Taq Gold (#4398833, Applied Biosystems, Monza, italy). After a first denaturation step at 94°C for 10 min, reactions were cycled 30 times (45 s at 94°C; 30 s at 58°C; 30 s at 72°C) with a final 10 min elongation step at 72°C.

Software and statistical analysis

All microsatellite data were first checked with the software micro-checker (Oosterhout et al., 2004) in order to spot null alleles and wrong size detections. The statistical analysis of the microsatellite data for the genetic variability measures, including the analysis of molecular variance (AMOVA) and the principal coordinate analysis (PCA), was performed using the latest version of the software Genalex 6.3 (Peakall and Smouse, 2006), while the Excel Microsatellite Toolkit (Park 2001) was used for calculating the polymorphism information content (PIC) for each allele.

A series of indices and parameters were calculated, such as allele frequencies and number of alleles, number of effective alleles, private alleles, expected and observed heterozygosis and fixation indices. The Fixation index (Wright’s inbreeding coefficient) was calculated using the software Genalex as (He - Ho)/He. Arcos-Burgos and Muenke (2002) gives an exhaustive review of this measure which, when applied to microsatellite data, describes the probability that a given locus becomes fixed. Populations showing high levels of homozygosity will have a Fixation index significantly different from zero.

All these parameters were calculated for the whole dataset and for the two populations: Suri and Huacaya. In order to assess if there were statistically significant differences between these results, the null hypothesis was tested by means of one way analysis of variance (ANOVA) for observed heterozygosis (Ho), expected heterozygosis (He), Fixation index (F), number of effective alleles (Ne) and PIC and by means of a non-parametric Mann- Whitney Test for number of alleles (Na).

The statistical population genetics package Fstat 2.9.3.2 (Goudet, 1995) was used to calculate deviation from Hardy-Weinberg equilibrium (HWE) and genotypic disequilibrium among loci applying a strict Bonferroni correction for multiple comparisons. The test for HWE and the test for linkage disequilibrium were carried out using 1300 randomizations: 1% table wide level was significant.

Cavalli-Sforza’s chord distance and Reynolds-Weir Cockerham distance were calculated using Gendist, an application of the software package Phylip version 3.69 (Felsenstein, 1989). The first measure is assuming a stepwise mutation model in an infinite allele model with equilibrium between mutation and genetic drift, whereas the second and the third measures are dimensional models assuming only genetic drift.

Finally, the genetic structure of the sample was investigated with the Structure 2.3.3 software (Pritchard, 2000) and the indentification of the most likely number of clusters (K) was made by the Evanno method (Evanno, 2005) using the online version of Structure Harvester (Earl, 2011). The burning period was set to 50,000 and repetitions of the MCMC chain to 106; the ancestry model chosen was the admixture model. Four replicates for each tested value of K (1–6) were performed.

Results

Hardy-Weinberg equilibrium and linkage disequilibrium

A number of indices have been calculated only on 12 microsatellites, excluding YWLL43, which is linked to the X sexual chromosome. Only one locus (LCA37) was found not to be in Hardy-Weinberg equilibrium after strict Bonferroni correction, showing excess of homozygosis (P<0.05). There was no difference in results when the sample was split into two different populations.

The test for genotypic disequilibrium between pairs of loci showed 9 loci associations out of 65 to be in some degree of linkage disequilibrium. The loci LCA8, LCA66 and LCA65 appear in 8 of the 9 associations showing linkage disequilibrium.

Heterozygosis, polymorphism information content and fixation index

When the dataset was considered as a single population, the average Ho for the 12 markers was high (Ho=0.766; SE=0.044), extremely close to the average He (He=0.769; SE=0.033), and the unbiased expected heterozygosis (UHe=0.778; SE=0.033), with an overall average fixation index of 0.004 (SE 0.036). The high mean number of alleles (Na=9.667; SE=0.77), effective alleles (Ne=4.89; SE=0.39) and the low fixation indices (F) confirm such high values of heterozygosis and genetic variability. The polymorphism information content (PIC) for each locus is in line with previous findings and ranges from 0.411 for locus YWLL46 to 0.826 for locus YWLL44. Table 2 shows Ho, He, UHe, PIC and F for each locus and as a mean for all loci with the relative standard errors.

Table 2 Population genetic parameters for a single population.

Locus	N	Na	Ne	Ho	He	UHe	F	PIC
YWLL46	45	6 (5)	1.75	0.44	0.43	0.43	-0.03	0.41
LCA65	44	13 (14)	5.50	0.90	0.81	0.82	-0.11	0.79
YWLL40	46	6 (7)	4.80	0.82	0.79	0.80	-0.04	0.75
LCA5	49	7 (13)	3.90	0.83	0.74	0.75	-0.12	0.70
LCA66	48	13 (24)	6.10	0.89	0.83	0.84	-0.07	0.81
LCA8	49	9 (14)	6.18	0.83	0.83	0.84	0.01	0.81
LCA99	48	11 (11)	4.03	0.64	0.75	0.76	0.14	0.72
YWLL44	49	11 (18)	6.41	0.77	0.84	0.85	0.08	0.82
LCA37	47	13 (19)	4.49	0.55	0.77	0.78	0.28	0.76
LCA94	47	7 (9)	4.09	0.68	0.75	0.76	0.10	0.72
YWLL36	49	10 (17)	6.16	0.91	0.83	0.84	-0.09	0.81
YWLL29	48	10 (9)	5.21	0.87	0.80	0.81	-0.08	0.78
YWLL43°	46	6 (10)	2.60	-	-	-	-	0.56
Mean	47.41	9.667	4.89	0.76	0.76	0.77	0.01	0.74
SE	0.48	0.77	0.39	0.04	0.03	0.03	0.03	0.03

N, number of individuals; Na, number of different alleles, in brackets values from previous studies (Obreque et al., 1999; Penedo et al., 1999; McPartlan et al., 1998); Ne, number of effective alleles; Ho, observed heterozygosity; He, expected heterozygosity; UHe, unbiased expected heterozygosity = [2N/(2N-1)] * He; F, fixation index; PIC, polymorphism information content;

° Ho, He, UHe and F not shown for X-linked YWLL43 locus.

There was no significant difference in results when the sample was analyzed as two separate populations (Huacaya and Suri). The null hypothesis was tested for Ho (P=0.69), He (P=0.61), F (P=0.95), Ne (P=0.69), PIC (p=0.61) and Na (P=0.37). Genetic variability parameters for the two populations are listed in Table 3. A number of private alleles were detected between the two phenotypes and a summary is given in Table 4 according to loci and phenotypes.

Table 3 Population genetic parameters for the two populations.

Population	Locus	N	Na	Ne	Ho	He	UHe	F
Huacaya	YWLL46	21	5	1.42	0.33	0.29	0.30	-0.12
	LCA65	21	9	4.34	0.95	0.77	0.78	-0.23
	YWLL40	21	6	4.47	0.85	0.77	0.79	-0.10
	LCA5	24	6	4.08	0.91	0.75	0.77	-0.21
	LCA66	23	11	5.68	0.91	0.82	0.84	-0.10
	LCA8	24	8	6.36	0.83	0.84	0.86	0.01
	LCA99	23	10	4.14	0.65	0.75	0.77	0.14
	YWLL44	24	9	5.78	0.70	0.82	0.84	0.14
	LCA37	24	9	4.15	0.50	0.76	0.77	0.34
	LCA94	24	6	3.61	0.62	0.72	0.73	0.13
	YWLL36	24	9	6.29	0.91	0.84	0.85	-0.09
	YWLL29	23	9	5.23	0.82	0.80	0.82	-0.02
	YWLL43§	22	5	2.28	-	-	-	-
Suri	YWLL46	24	5	2.11	0.54	0.52	0.53	-0.02
	LCA65	23	11	5.71	0.87	0.82	0.84	-0.05
	YWLL40	25	6	4.92	0.80	0.79	0.81	-0.01
	LCA5	25	6	3.37	0.76	0.70	0.71	-0.08
	LCA66	25	9	6.28	0.88	0.84	0.85	-0.04
	LCA8	25	9	5.68	0.84	0.82	0.84	-0.01
	LCA99	25	8	3.85	0.64	0.74	0.75	0.13
	YWLL44	25	10	6.06	0.84	0.83	0.85	-0.01
	LCA37	23	13	4.76	0.60	0.79	0.80	0.23
	LCA94	23	7	4.40	0.73	0.77	0.79	0.04
	YWLL36	25	8	5.98	0.92	0.83	0.85	-0.10
	YWLL29	25	10	5.04	0.92	0.80	0.81	-0.14
	YWLL43°	24	5	2.75	-	-	-	-
	N	Na	Ne	Ho	He	UHe	F
Huacaya	Mean	23	8.08	4.63	0.75	0.74	0.76	-0.01
	SE	0.36	0.54	0.39	0.05	0.04	0.04	0.05
Suri	Mean	24.41	8.50	4.85	0.78	0.77	0.79	-0.01
	SE	0.26	0.66	0.36	0.03	0.02	0.02	0.03

N, number of individuals; Na, number of different alleles; Ne, number of effective alleles; Ho, observed heterozygosity; He, expected heterozygosity; UHe, unbiased expected heterozygosity; F, fixation index;

° Ho, He, UHe and F not shown for X-linked YWLL43 locus.

Table 4 List of private alleles in each population.

Population	Locus	Alleles
Huacaya	YWLL46	95
	LCA65	175, 179
	LCA5	190
	LCA66	220, 244, 252, 256
	LCA99	278, 290, 294
	YWLL44	108
	YWLL43	148
	YWLL36	168, 172
Suri	YWLL46	111
	LCA65	165, 177, 185, 189
	LCA5	192
	LCA66	240, 260
	LCA8	245
	LCA99	268
	YWLL44	116, 134
	YWLL43	162
	LCA37	130, 138, 148, 170
	LCA94	197
	YWLL36	170
	YWLL29	238

Genetic distance, AMOVA, PCA and genetic structure

Genetic distance calculated by Pairwise Population Matrix of Nei’s Genetic Distance and Unbiased Nei’s Genetic Distance (Nei, 1978) were 0.062 and <0.0001, respectively. When calculated as the Cavalli-Sforza’s chord distance (Cavalli-Sforza and Edwards, 1967) results showed 0.03. Reynolds-Weir Cockerham distance (Reynolds et al., 1983), suggested to be more precise in the calculation of genetic distances between closely related species and breeds (Laval et al., 2002), was 0.04.

The PCA calculated on the distance matrix among the individuals in FIgure 1 graphically shows how samples from the Suri and Huacaya datasets overlap and do not segregate into different groups. The first 3 dimensions of the PCA explain 64.64% of the total variance. All the variance observed in the two populations with the AMOVA test, calculated both by Fst and Rst values, was due to variation within populations (100%) and not to variation between populations (0%). The software Structure rendered a maximum likelihood for K=2 both with the classic method by Pritchard (2000) and with that by Evanno (2005), therefore suggesting a possible dual ancestry for the analyzed sample. However, individuals from both phenotypes were assigned in equal proportions to the two clusters without a clear-cut distinction between the two groups, as shown in FIgure 2.

Figure 1 Principal coordinates analysis of Suri and Huacaya based on the distance matrix from molecular data. Red, Huacaya; green, Suri.

Figure 2 Bar plot for cluster assignment. 1, Huacaya; 2, Suri. Colored bars indicate the probability of assignment to either cluster 1 (green) or cluster 2 (red).

Discussion

In terms of genetic variability within the Peruvian alpaca sample analyzed, all parameters and findings, such as number of alleles (Na=9.667; SE=0.772), number of effective alleles (Ne=4.89; SE=0.388), observed and expected heterozygosis (Ho=0.766, SE=0.044; He=0.769, SE= 0.033), show that the Peruvian population is still conserving high genetic variability and does not show any sign of artificial selection pressure for the studied loci. The low fixation indices for these loci confirm this interpretation of the data (F=0.004; SE=0.036) and suggest that the microsatellite panel used is suitable for genetic diversity studies. In comparison with the data collected and analyzed by the ARI since 1998 and by several other authors, the sample in this study shows the presence of a high number of alleles, matching the whole allelic range described in previous publications (Obreque et al., 1999; Penedo et al., 1999; McPartlan et al., 1998). Values for the polymorphic information content (PIC=0.746, SE=0.033) were also in line with previous findings. Only one locus (YWLL46) showed a PIC value less than 0.7 (PIC=0.411), which reflects the lower than average Na and Ne found for this specific locus (Na=6; Ne=1.753) in previous studies. Nevertheless, this value is higher than that from previous bibliographic data (Lang et al., 1996).

There was no significant difference between the two phenotypic groups when these parameters were evaluated separately. In terms of genetic distance and differentiation between the two phenotypes the PCA analysis shows no separate segregation or grouping of Suri and Huacaya individuals.

To our knowledge, no previous studies have investigated genetic distance between Suri and Huacaya alpaca by means of codominant markers, although microsatellite markers have been used in the past to investigate the phylogeny of a wider group of South American Camelids (Wheeler et al., 2006, Bustamante et al., 2002).

When ANOVA was carried out considering the two phenotypes as two separate populations, it clearly identified the source of all variance in the component within populations, excluding any source of variance to be found between populations. This result is supported by Nei’s index of genetic distance, Cavalli-Sforza’s chord distance and Reynolds-Weir Cockerham distance, which also show no differentiation between the two populations. There are two important factors to be taken into consideration while interpreting these results. Firstly, the germplasm established 20 years ago at the experimental station of Quimsachata is not subjected to genetic selection and was created for the sole purpose of conserving the genetic diversity of the species. Secondly, although Suri and Huacaya alpacas at the experimental station have been bred and managed separately since the creation of the germplasm, the time interval of 20 years is unlikely to generate genetic differentiation between the two phenotypes, especially considering the absence of selection, the reproductive physiology of the species and its generation time (Mason, 1973).

Conclusions

Given these considerations, if a secondary breed structure had been present within the species at the time of the creation of the germplasm, it would have been preserved by the breeding practices in place at the experimental station. Nevertheless, the data obtained from the 13 loci suggest no genetic divergence between the two phenotypes and do not support the idea of two distinct populations of Peruvian Suri and Huacaya alpacas. Furthermore, the two phenotypes have similar genetic parameters in terms of allelic frequencies and genetic variability, showing high values both in terms of allelic richness and heterozygosis.

Acknowledgment:

the authors would like to thank the LGS laboratories for carrying out the multiplex amplifications and band scoring.