Introduction

The alpaca is one of the two domestic South American Camelids and is raised primarily for fiber production. Alpaca fiber is composed of alpha keratins and keratin-associated proteins (KAPs) (^{Powell & Rogers 1997}). The latter are a complex group of proteins, which have high levels of cysteine, or glycine and tyrosine. It is hypothesized that these proteins play an important role in defining the physical-mechanical properties of the fiber, therefore sequence variation and variation in the expression of the genes encoding KAPs may be associated to variations in the quality of the fiber (^{Gong et al. 2012}, Powell & Rogers 1997). KAPs are divided into two main groups, those with high sulfur content and those with high glycine and tyrosine content. The first group can be separated into two subgroups, high-sulfur proteins with less than 30% cysteine residues and very high-sulfur proteins with more than 30% cysteine residues (^{Gillespie & Broad 1972}, Gillespie 1972).

The genetic improvement of alpacas through selective breeding is a subject of interest for breeders and the textile industry involved in alpaca fiber processing. Breeders aim to augment their revenue streams by augmenting both the quantity of fiber produced and diminishing the diameter of individual fleece fibers per sheared animal. Meanwhile, the textile industry endeavors to procure a larger supply of extra-fine and fine-quality fiber, which will subsequently be utilized in the fabrication of textiles with higher economic value in the market (^{Gutiérrez 2008}).

Single nucleotide polymorphisms (SNPs) are characterized by single-base variances in the DNA sequence among individuals within a species. These molecular markers serve as indispensable tools in the genetic improvement of various animal species through genomic selection. SNPs hold significant value in association studies owing to their extensive variability and widespread distribution across the genome. Additionally, they play a crucial role in establishing kinship relationships, as the majority of SNPs exhibit uniqueness within a given population (^{Pierce 2006}).

From previous studies on SNPs in alpacas, it is important to mention that according to ^{Jones et al. (2019}), the gray color in alpacas is the result of a SNP c.376G>A (p.Gly126Arg) located in exon 3 of the KIT gene. ^{Likewise, Salas (2019}) found 34 SNP-type variants in genes involved in the structure and development of alpaca fiber (KRT31, HOXC13 and EDAR); however, no association was found between the gene fragments analyzed and the diameter of the alpaca fiber. ^{Foppiano (2016}) identified 27 SNP markers in five genes of the keratin-associated protein family (KRTAP1-2, KRTAP6-1, KRTAP9-2, KRTAP11-1, KRTAP13-1), and conducted an association study between fiber diameter and two SNPs of the KRTAP11-1 gene but did not find a significant association in a sample of 152 white Huacaya alpacas from Puno (Peru). ^{Fernandez et al. (2019}) performing a bioinformatics analysis based on orthologous genes from sheep, goat, and human, identified 48 SNPs in 22 alpha keratin genes in alpacas. Finally, ^{Calderon et al. (2021}) reported the development of a 76k alpaca SNP microarray, which included 302 SNPs that were found in candidate genes for fiber quality and color.

The objective of this study was to identify single nucleotide polymorphism molecular markers within genes encoding keratin-associated proteins (KRTAPs) in alpacas. Quality parameters such as the absence of SNPs in the flanking sequence (70 bp), minor allele frequency, rate of samples with analyzed sequences, and Illumina Score were considered to identify reliable SNPs suitable for inclusion in an alpaca SNP microarray.

Material and methods

Genomic data. For sequence comparisons to identify SNPs we used the alpaca reference genome (VicPac3.1, GenBank accession number: GCA_000164845.4) and nine alpaca genomes available in NCBI. Three of these were generated by the Universidad Peruana Cayetano Heredia (Accession number: PRJNA340289) and six genomes by the Universidad Nacional Agraria La Molina (Accession number: PRJNA685331) (Anex 1). We also used sequence reads from 150 reduced representation DNA libraries with DNA fragments generated after digestion with the ApeKI enzyme and 150 generated after double digestion with the PstI and Msp enzymes (^{Calderon et al. 2021}).

The sequences, number of exons and location (scaffold and chromosome) of 34 keratins associated proteins genes reported in alpacas were obtained from the NCBI. To determine the initial and final position of each exon, as well as of the introns and untranslated regions, the Sequence Viewer (https://www.ncbi.nlm.nih.gov/projects/sviewer/) was used.

Identification of single nucleotide polymorphisms. Having obtained the positions of each exon, a coordinate file (bed file) of the exons, introns, and untranslated regions (UTRs) of each gene was generated. To align the sequencing data of the nine alpaca genomes, referred to above, to the 300 alpaca reduced representation DNA libraries, also described above, and to the VicPac3.1 reference genome we used BWA v.0.7.10 (^{Li & Durbin 2010}). While single nucleotide polymorphisms (SNPs) were identified using BCFtools (https://samtools.github.io/bcftools/bcftools.html),

Selection of SNPs. A reliable SNP should be devoid of other SNPs in its flanking sequences (35 bp 5’ of the SNP and 35 bp 3’ of the SNP, for a total of 71bp). Selection of SNPs for inclusion in a microarray requires a minor allele frequency (MAF) of ≥ 0.05 with a genotyping rate > 45 %, and an Illumina Score ≥ 0.6 (^{Laurie et al. 2010}). The absence of SNPs within flanking regions was assessed by creating a vcf file with the genotype information using the KGD v0.8.2 software with default parameters (https://github.com/AgResearch/KGD). MAFs and genotyping rate for each SNP was calculated using the PLINK software and the Illumina Score was calculated using the Illumina Design Studio (https://designstudio-array.illumina.com). Thus 35 KRTAPs SNPs were included in the 76K alpaca SNP microarray (^{Calderon et al. 2021}).

Validation of SNPs. The selected 35 KRTAPs SNPs were validated by genotyping 937 white Huacaya alpacas (58 males and 879 females) with the 76K alpaca SNP microarray. Genotyping was done at NEOGEN GeneSeek Laboratories (Nebraska, United States). For this study, only the genotyping results of SNPs located in KRTAPs were analyzed. Quality control was performed using the PLINK v1.90p program (^{Chang et al. 2015}). Samples with a genotyping rate ≥ 90% and SNPs with a genotyping rate ≥ 90% were retained.

Results

The location and sequences of the 34 KRTAP genes annotated in the reference genome of alpacas were obtained, of which according to the classification of ^{Gong et al. (2012}), 19 are KRTAPs encoding high sulfur content proteins, 6 are KRTAPs encoding ultra-high sulfur content proteins, and 9 are KRTAPs encoding high glycine and tyrosine content proteins (Table 1).

Table 1 KRTAPs annotated in the reference VicPac 3.1 alpaca genome. 

A total of 67 SNPs were identified (Table 2), of which 48 are distributed in 18 KRTAPs encoding high sulfur content proteins, 8 in 5 KRTAPs encoding ultra-high sulfur content proteins, and 11 in 9 KRTAPs encoding high glycine and tyrosine content proteins.

Table 3 shows the results of the validation of 35 SNPs in 937 alpacas that were genotyped. After quality control, 936 alpacas and all SNPs were retained. Three SNPs were monomorphic (MAF=0), nine had MAFs between 0.006 and 0.049, and 23 had MAFs ≥ 0.05.

Discussion

Sixty-seven SNPs were identified for the 34 KRTAPs genes that are annotated in the alpaca reference genome. Of these, 35 SNPs are included in the 76K alpaca SNP microarray (^{Calderon et al. 2021}) and 32 were confirmed as true SNPs and three were monomorphic in the alpaca population that was analyzed. These 32 validated SNPs are localized in 14 KRTAP genes. The three monomorphic SNPs indicate that the polymorphisms identified are not present in the alpaca population used for our analysis and/or are of low frequency in general. The latter can be true for any SNP because their frequencies will vary depending upon the sample populations being analyzed unless their frequencies are closer to 50%.

The constraints of this study lie in the potential for gene locations to shift with advancements in the alignment of the reference genome utilized. Furthermore, enhancing genome assembly is imperative for accurately determining the physical locations of genes and, consequently, single nucleotide polymorphisms (SNPs). Additionally, the presence of sequencing errors in the genomes and reduced libraries are plausible. Utilizing libraries is crucial as it ensures that sequenced fragments possess adequate repetitions, thereby facilitating the generation of a more robust consensus sequence.

In the KRTAP11-1 gene, a single nucleotide polymorphism (SNP) (g.12376109T>G) was detected, consistent with SNP33 as identified by ^{Foppiano (2016}). Foppiano (2016) additionally documented five additional SNPs within the same gene, though these variants showed no association with fiber diameter. Our study, however, did not detect all of the SNPs reported by Foppiano (2016), possibly due to our stringent filtering criteria based on quality parameters outlined earlier. Foppiano (2016) employed only allele frequency as a filtering criterion.

The KRTAP13-1 gene consists of three copies (LOC102533350, LOC102532866, LOC102532037). The SNP (g.12846534G>C) identified within KRTAP13-1 (LOC102532866) matches the one reported by ^{Foppiano (2016}), who noted its conversion from Tryptophan to Cysteine in the protein sequence, a finding consistent with our study's results. Notably, Foppiano (2016) aligned the alpaca genome (Vicugna-Pacos 2.0.1) with the human KRTAP13-1 gene sequence to ascertain the identified alpaca gene sequence. Moreover, Foppiano (2016) identified three additional SNPs within the same LOC102532866 sequence that our study did not detect due to stringent quality SNP identification filters. However, our study successfully identified SNPs within all three copies of the KRTAP13-1 gene as annotated in the NCBI. We identified four SNPs in LOC102533350, two of which are non-synonymous, resulting in alterations in the protein sequence from Threonine to Serine (g.12923905C>G) and Valine to Phenylalanine (g.12924045G>T). Similarly, we observed two SNPs in LOC102532037 leading to changes in the protein sequence, from Alanine to Serine (g.12859532G>T) and Serine to Cysteine (g.12859761C>G). In goats, ^{Fang et al. (2010}) identified a SNP (T>G) within the gene's coding region, significantly associated with fiber diameter in Xinjiang (cashmere) breed goats. In our study, we identified a G>T SNP at positions 12924045 and 12859532 of the same gene. We recommend employing these SNPs within KRTAP13-1 (LOC102533350) for association analyses with fiber characteristics in alpacas.

In KRTAP24-1 gene, four SNPs were identified and validated, three of which are non-synonymous SNPs. One produces a change from Alanine to Valine (g.12980054C>T), which is similar to the SNP (c.656T/C) identified by ^{Wang et al. (2019}) in goats. This SNP causes a change of Valine to Alanine in the protein sequence and was reported to affect the fiber diameter in cashmere goats, while the other two changed Glycine to Histidine (g.12979659A>T) and Asparagine to Serine (g.12979802A>G) were not reported in previous studies in other species. The fourth SNP identified in this study is synonymous (g. 12979560 T>G). However, it would be important to consider all of them for an association study analysis on fiber characteristics in alpacas.

For the KRTAP26-1 gene, four SNPs were identified, one synonymous that was validated and three non-synonymous, not validated, that produce the following changes in the protein sequences, Serine to Proline (g.12947455T>C), Leucine to Glutamine (g.12947708T>A), and Arginine to Cysteine (g.12947881C>T), respectively. Similarly, ^{Li et al. (2017}) identified four SNPs in this gene but the one that was associated with lower mean wool diameter, standard deviation of wool diameter, mean comfort factor, and higher wash performance in Merino x Southdown sheep (c.277A/G) changes Serine to Glycine, was not found in our study. Again, it would be important to analyze the association with fiber characteristics of these SNPs in alpacas.

A single nucleotide polymorphism (SNP) within the untranslated region (UTR) of the KRTAP9-2 gene was identified. Additionally, ^{Foppiano (2016}) reported a SNP in the 3'UTR; however, the former SNP was excluded from further analysis due to its failure to meet the allelic frequency criterion (≥ 5%), while the latter SNP was not detected in our study. Foppiano (2016) also noted two additional SNPs within the exonic region, which were not included in our analysis due to their failure to meet the quality parameters required for inclusion in the 76K SNP microarray. In studies involving cashmere-producing goats, ^{Yu et al. (2008}) identified a 30-nucleotide insertion/deletion variation in the exon of the KRTAP9-2 gene, while ^{Wang et al. (2012}) identified a C/T SNP variation in the same gene. Furthermore, in a separate investigation, Wang et al. (2014) conducted a gene expression profiling study across different phases of the fiber cycle, revealing significant differences in the expression levels of this gene between goats exhibiting high and low production of cashmere fiber. Consequently, further investigation into the presence of additional molecular markers, such as insertion/deletion mutations, remains warranted in alpacas.

In the single exon of the KRTAP7-1 gene, we found a non-synonymous SNP that changes a Serine to Arginine (g.12420218C>G), comparatively, ^{Liu et al. (2014}) identified a non-synonymous SNP in Chinese merino sheep, but this changed Asparagine to Lysine in the protein sequence. In Liaoning goats, ^{Jin et al. (2011}) found that the expression of this gene was 1.28 times higher in secondary follicles than in primary follicles, suggesting that it probably plays an important role in regulating fiber diameter.

For the KRTAP8-1 gene, we found two SNPs, one synonymous (g.12436062T>C) and the other non-synonymous (g.12436103G>A), unlike ^{Liu et al. (2011}), who found two SNPs in goats, but both were synonymous mutations, and one (T113G) was significantly associated with cashmere fleece weight, but not with fiber diameter. Our non-synonymous SNP is included in the 76K alpaca SNP chip and would be important to consider it for association studies with fiber characteristics.

For the KRTAP20-2 gene, we identified five SNPs, one of which is non-synonymous and changes an Alanine to Valine (g.12678801C>T), while the other four identified SNPs are synonymous. ^{Bai et al. (2018}) also reported a non-synonymous SNP (c.160A>T), it generates a premature terminal codon associated with mean diameter and wool curvature in sheep (Merino x Southdown).

In the single exon of the KRTAP21-1 gene we did not find any exonic SNPs, while ^{Li et al. (2019}) reported four exonic SNPs in Merino x Southdown lambs and found one associated with a positive variation in wool washing performance. However, two SNPs were identified in the UTR of the KRTAP21-1 gene, which might be involved in gene expression regulation.

The validation study of 35 out of the 67 identified SNPs in KRTAP genes was conducted successfully, as they were incorporated into the 76K alpaca SNP microarray. Alternative approaches will be employed to validate the remaining SNPs, such as gene sequence comparisons among genomes as they become accessible (^{Foppiano 2016}, Palloti et al. 2023).

To summarize, 32 out of the 67 identified SNPs have been validated in a population of 936 white Huacaya alpacas. While the incorporation of these SNPs into the 76K alpaca SNP chip will aid in elucidating their association with fiber characteristics, it is equally imperative to explore the biochemical attributes of these proteins and their follicular expression levels in alpacas in future whole-genome association studies.