A Novel SNPs of KISS1 Gene Strongly Associated with Litter Size in Indonesian Goat Breeds

Kisspeptin is a protein encoded by the KISS1 gene which behaves as a key role by stimulating gonadotropin releasing hormone (GnRH) neuron activity directly in the reproductive axis. The objective of the present study was to determine the genetic diversity within intron 1 of KISS1 gene and to verify their association with fecundity traits which can be devoted as a marker assisted selection (MAS) for breeding selection. This study was established on three Indonesian native goat populations (Kacang, Kejobong, and Senduro). The PCR products were then sequenced in both directions. The DNA sequencing alignment resulted fifteen variants (one indel and fourteen SNPs), with SNP1, SNP2, SNP3, and SNP9 being novel single nucleotide polymorphisms (SNPs) of KISS1 gene intron 1. Genotype, haplotype, and parity were significantly associated with litter size. SNP8, SNP9, and SNP10 were correlated with overall means of litter size (LS) and LS at the first and the third parities (p<0.05). Additionally, novel SNP9 were in remarkably strong linkage disequilibrium with SNP8 and SNP10 (D’=1.00; r≥0.58; χ 2 ≥13.38), and allele A had putative binding sites for the NFIC which plays an important role in activating the expression of KISS1 gene. The H2 haplotype (CATAGCGCAACGCT) was found to have the highest litter size (p<0.0001). CC genotype at SNP8, AA genotype at novel SNP9, GA genotype at SNP10, and H2 haplotype were the excellent genotypes and haplotype that associated with the superior LS (p<0.05). Therefore, this result led to presume that these three SNPs and H2 haplotype can be acknowledged as prominent genetic markers for goat-breeding selection.


INTRODUCTION
Goats were reared for several reasons, including their abilities to adapt in diverse ecosystems (frigid, harsh, and arid zones), adaptive under low input environments, tolerance to heat stresses, diseases, and parasites, and their abilities to provide socioeconomics, cultures, and religious benefits, as well as cash readily for immediate needs and meaningful income among poor households.Indigenous goats have been part of rustic subsistence in many developing countries (Miller et al., 2012;Kaumbata et al., 2020).Nonetheless, native goats have been extensively reared in Indonesia, resulting in a negative genetic correlation among adaptation and production traits.Therefore, it might be hard to compromise due to selection for production efficiency (Gama & Bressana, 2011).
The goat population in Indonesia is nearly nineteen million, contributing 72.553 tons of meat per year from 7.144.010heads (FAOSTAT, 2019).The Indonesian Ministry of Agriculture (2020) reported that goat meat redounds 1.8% of national meat productions.The import of goat meat products reached 1.49%, and total goat imports were 387.663 heads.This condition demonstrated that goat population's scale-up remains critical to meet the national needs.
Domestic animal genetic diversity is essential for meeting current production needs in various environments in the future and changing breeding objectives, and long-term genetic improvements (Liu et al., 2019) through enhancement of the reproduction rate of goat breeds (Ahlawat et al., 2015).Nevertheless, the progress in reproductive ability traits is low due to the long generation interval (Rupp et al., 2016).On the other hand, native breeds/populations are critically threatened to become extinct without recorded phenotypic and genotypic data (Periasamy et al., 2017).Traditional breeding programs and/or selection procedures using genetic information obtained through marker-assisted selection could improve the prolific traits, economic gain, and genetic gain (El-Tarabany et al., 2017).
Identifying which genes are responsible for high fecundity traits/litter size (LS) is crucial.According to genetic research, KISS1 gene is a major fecundity gene in goats.Polymorphisms of KISS1 gene were associated with higher LS (Cao et al., 2010;An et al., 2013).In mammals, kisspeptins (Kp) are encoded by the KISS1 gene and play an important role in reproductive traits.Kp has been identified as pivotal aspects and crucial gatekeepers of reproductive action and maturation, puberty onset, sexual differentiation of the brain, adult regulation of gonadotropin (GnRH) secretion, and the metabolic control of fertility (Pinilla et al., 2012).KISS1 gene in goats has been located in chromosome 16.KISS1 gene consists of two coding regions (exons) and one single non-coding region (intron), and the transcript length is 408 bp and encodes 135 amino acids.This gene reaches around 2.62 kilo bases (ENSCHIT00000037363.1).
Previous research in Indonesia characterized the polymorphisms of KISS1 gene intron 1 and exon 1 of Ettawa Grade (EG), although the researchers could not find the association with LS and other reproductive traits (Mulyono et al., 2019;Hardyta et al., 2020).However, to make the best use of limited selection funds, specific indigenous breeds must be prioritized for selection.Therefore, further research in other fragments of KISS1 gene is needed in Indonesian native goat breeds.Three goat breeds in Indonesia that are wellknown and registered at the nationwide level, namely Kacang (KC) and Kejobong (KJ) goats, are meat-purpose breeds, while Senduro (SD) goat is a dual-purpose breed.The KC goat is the indigenous goat in Indonesia and Malaysia, commonly known throughout South East Asia.They have a high fecundity and heat tolerance in harsh environments despite limited potential growth (Tsukahara et al., 2008).KJ goat is a crossbreed between KC and EG goats (Lestari et al., 2018), whereas SD was identic with EG goats (Ciptadi et al., 2019).
No studies have found a strong linkage between SNPs associated with LS of KISS1 gene in Indonesian native goat breeds.Hence, this research was conducted to determine the genetic variations of intron 1 of KISS1 gene in Indonesian native goat breeds and examine the relationship of its SNPs and haplotype with fecundity traits that could be used as a marker-assisted selection (MAS) for the genetic breeding selection.

Selection of Animals, Collection of Data, and Biological Substances
In the current research, 90 female goats (does) were selected based on the breed/phenotypic characteristic recognized from the specific breed, doe's parity (the first to the fifth parities), LS, and farms from different regions.The samples were collected from February to June 2020.Three different breeds of goats were used, i.e., Kacang (KC), Kejobong (KJ), and Senduro (SD), with 30 samples each.The other research on gene polymorphism associated with economic traits on goats used fewer samples (Abdel-Aziem et al., 2018;Antonius et al., 2020;Sasi et al., 2020;Rahmatalla et al., 2020;Rossi et al., 2017).They were reared by three groups of farmers in their native environment and uncorrelated.All the animals used at the present study were under the standard rule of animal treatment designated by the Animal Ethics Committee of the Faculty of Animal and Agricultural Sciences, Diponegoro University (No. 57-06/A-4/KEP-FPP).The farmers implicated in the study were notified and approved to use the animal in the present research.
Five milliliters of whole blood samples were collected from the jugular vein and stored in a sterile EDTA vacutainer tube at -20 °C.The genomic DNA was extracted using the GeneJET Genomic DNA Purification Kit (Thermo Scientific, USA), conforming to the manufacturer's guidance.The quantity and quality of genomic DNA extraction were established using agarose gel electrophoresis (1%) and Nanodrop spectrophotometer Uvidoc HD6 (UVItec Ltd., Cambridge, UK).The electrophoresis showed a clear single band on agarose 1%, and the OD 260/280 ratio ranged ≥1.7 (according to the kit protocol), demonstrating a good quality of DNA extraction.
PCR cycling program consisted of initial denaturation at 95 °C for 5 min, followed by 35 cycles of denaturation at 94 °C for 30 s, annealing at 60 °C for 30 s, extension at 72 °C for 30 s, and final extension at 72 °C for 7 min.Amplification products were defined by electrophoresis in 2% agarose gel (Invitrogen, Life Technologies Co, CA) at 100V for 30 min alongside 100 bp DNA marker (Geneaid, 1 st BASE).The bands were captured under UV light and the gel appearance was documented using digital gel imager equipment (UV Transilluminator, Uvitec Cambridge).The KISS1 gene-PCR developed 1061 bp of DNA fragment, as shown in Figure 1.The PCR products (n=90) were sequenced in two directions using commercial service.Due to the poor DNA sequences, only 80 nucleotide sequences were used to screen the variation of nucleotides.

Statistical and Data Analysis
SNP genotyping.The multiple nucleotide sequences were aligned by software MEGA X 10.1 version (Kumar et al., 2018) and compared with the sequence from the GenBank assembly accession No: GU142847.1 using BLAST program (www.ncbi.nlm.nih.gov/blast/).The recognized SNPs were then combined to establish haplotypes (Swain et al., 2020).
Linkage disequilibrium and association analysis.The number of haplotypes, linkage disequilibrium (D'), and correlation coefficient (r) test were calculated by DnaSP 6 th version (Rozas et al., 2017).The significance of ꭓ2 value was compared using the Bonferroni correction for multiple tests within D' (Fu et al., 2014).Genotype and allele frequencies were calculated using the direct counting methods following Nei & Kumar (2000): where x ii was the genotype frequency of genotype ii, n ii was the number of samples with genotype ii, and n was the number of total samples.
where x i was the allele frequency of allele i, n ii was the number of samples with genotype ii, n ij was the number of samples with genotype ij, and n was the number of total samples.Expected heterozygosity and Hardy Weinberg Equilibrium (HWE) were estimated using Arlequin Software Package version 3.5 (Excoffier et al., 2010).Association between individual SNPs and haplotypes with reproductive traits in three goat breeds were investigated using fixed General Linear Model (GLM) of SAS Software (SAS ® University Edition, 2018).Multiple comparisons of the means were analyzed using Tukey's post hoc.
The fixed model used for prolific traits is: where y ijklm is the observed value of a dependent variable (LS); µ is the overall mean of the population; G i is the fixed effect of ith genotype (i = 1, 2, 3); H j is the fixed effect of jth haplotypes (j = 1, 2,…,14); F k is the fixed effect of kth farm group (j = 1, 2, 3); P l is the fixed effect of lth parities (l = 1, 2, 3, 4, 5); B m is the fixed effect of mth breeds (m = 1, 2, 3) and e ijkl is the random residual for y ijklm .When there was p<0.05, it was examined significant statistically.

PCR and Analysis of Intron 1 KISS1 Gene Sequence
According to the electrophoresis results, a single fragment of 1061 bp in length corresponded to the main part of KISS 1 gene at intron 1.The fragments were separated using 2% agarose gel (Figure 1), indicating that amplification product size was consistent with the target and had a good specificity, allowing for immediate sequencing in both directions.
The sequence alignment analysis revealed one insertion/deletion and fourteen polymorphic sites in the intron 1 KISS1 gene on three Indonesian native goat breeds (Table 1).The sequences of different genotypes are shown in the appendix.The nucleotide's numbering was based on the goat KISS1 gene sequence (GenBank accession number: GU142847.1)related to intron 1.
It is proper to address that one indel and five SNPs (SNP4, SNP6, SNP11, SNP12, and SNP13) out of fifteen mutations recognized in recent studies conform to prior research (Cao et al., 2010;An et al., 2013).The current study detected four novel SNPs of KISS1 gene at intron 1 in Indonesian native goat breeds, which are SNP1, SNP2, SNP3, and SNP9.

Allele Frequencies, Genotype Frequencies, Expected Heterozygosity (He), and Hardy Weinberg
Equilibrium (HWE) Based on the nucleotide sequence of KISS1 gene intron 1 fragments, the allele and genotype frequencies of each SNPs in 3 Indonesian native goat breeds were defined.All the three possible genotypes were presented in most of the four novel SNPs in KISS1 gene.Allele and genotype frequencies were given in Table 2.The ten mutations (g.1960-1977 DelTATGAGTAGCCCCCCTGC, SNP1, SNP7, SNP8, SNP9, SNP10, SNP11, SNP12, SNP13, and SNP14) were genotyped in three goat breeds.However, polymorphisms at SNP1, SNP8, and SNP12 could not be detected in SD goat (monomorphic).Polymorphisms at SNP13 were discovered only in KJ goats.Therefore, polymorphism at SNP14 and wild type (II) at g.1960-1977 DelTATGAGTAGCCCCCCTGC was found in SD goats only.
Hardy Weinberg Equilibrium (HWE) and expected heterozygosity (He) were calculated based on the genotype frequencies in all breeds (Table 3).Several deviations in HWE were found between breeds in the same site.
In all populations, we found that genotypes at one indel and four SNPs (g.19601977DelTATGAGTAGCCCC CCTGC, SNP8, SNP9, SNP11, and SNP12) were not fit

Association Analysis Based on Genotype
The overall mean of LS was calculated 2.23±0.09with a range from one to five.The mean LS for KC, KJ, and SD are 2.17±0.65,2.33±0.84,and 2.2±0.96,respectively, and the frequency of the birth type of each goat breed was given in Table 4.The overall frequency of single kids was 18.89%, twins was 43.3%, triplets was 34.44%, and quadruplets or more was 3.33%.
The correlation study among LS and fixed effects (goat breeds, parity, haplotypes, farms, and genotypes) were done using fifteen variants of KISS1 gene at intron 1 in three goat breeds.In the current research, the goat breed and farm have no significant correlation with LS (Table 5).
The parity had a significant association with LS (p<0.0001).The greatest LS was found in the third to the fifth parities, while the lowest LS is found in the second parities (Table 5).Furthermore, due to the lack of sample in the fourth and the fifth parities, further statistical analysis was conducted using the first to the third parities.The least square means with standard error of litter size for different genotypes and parities are shown in Table 6.The overall mean for the first parity, the second parity, and the third parity were 2.67±0.45,1.88±0.09,and 3.00±0.00,respectively.LS at the third parity differed significantly from LS at the second parity.CC genotype at SNP1 and SNP13 correlated significantly with higher LS at the first and second parities (p<0.05).However, there was no significant association between genotypes at g. [1960][1961][1962][1963][1964][1965][1966][1967][1968][1969][1970][1971][1972][1973][1974][1975][1976][1977] DelTATGAGTAGCCCCCCTGC and SNP14 at the first and the second parities.
The association analysis between different genotypes and LS on three Indonesian native goat breeds (Table 7) showed no significant difference in five SNPs (SNP2, SNP3, SNP4, SNP5, and SNP6),

Association Analysis of Haplotypes
Strong linkage disequilibrium (D') between SNPs of KISS1 gene intron 1 has been identified in Indonesian native goat breeds.The D' value range from 0.042 to 1, and coefficient of correlation (r) range between 0.01 to 1 (Table 8).The DnaSP software distinguishes 14 haplotypes in intron 1 of the KISS1 gene using 15 SNPs in the non-coding region from three goat breeds.The haplotypes were discovered to correlate significantly with LS (p<0.0001).The goats having haplotype 2 (CCATAGCGCAACGT) had the highest LS (4±0.33), while haplotype 11 (CATTGCACAGTGCT) had the lowest LS (1±0.31).Association among haplotype polymorphisms of KISS1 gene with LS was available in Table 9.

TFBS Prediction of SNP9
In the current study, the prominent novel variant (SNP9) correlated with higher litter size was used to predict transcription factor binding sites for the intron region of KISS1 gene of goats.The VEP analysis predicted that this mutation has a transcript feature type and a protein-coding biotype, indicating a regulatory role in protein expression.Animal TFDB was used to screen any TFBS alteration due to the native and mutant variants of the investigated fragment of the KISS1 gene.In the intron 1 of the KISS1 gene sequences where the A/G polymorphism was discovered, two transcription factors, NFIC and ESR1, were discovered to bind (Table 10).
The previous study about the polymorphisms of KISS1 gene as a candidate gene for prolific traits has been explored in goats.In Jining Grey goat, there are three mutations in intron 1 (G296C, G454T, and T505A), 18 bp insertion/deletion (indel) variant in intron 1 and no mutation found in exon 2, further two mutations (G3433A and C3688A) were located in exon 3 (Cao et al., 2010).Another study could detect ten polymorphisms in KISS1 gene in three Chinese goat breeds (g.384G>A, g.1147T>C, g.1417G>A, g.1428_1429delG, g.2124C>T, SNP6, SNP11, SNP12, SNP13, g.3864_3865delCA, and g.3885_3886insACCCC (An et al., 2013).TT genotype at intron 1 had a higher progesterone level in comparison with TA genotype in Damascus and Zarabi goat breed (El-Tarabany et al., 2017).Ferreira et al. (2020) argued that the SNPs that affected a production trait could be used as MAS for selection programs in small herds or when genomic selection is not possible or adjustable.
Another study in the coding region on KISS1 gene gave a contradictory result.According to Feng et al. (2009), exon 1 of the KISS1 gene had no polymorphisms.Similar to previous research, Mulyono et al. (2019) could not discover polymorphisms in exon 1 of Etawah Grade goat.On the contrary, Mekuriaw et al. (2017) found a polymorphism in exon 1 and exon 2 KISS1 gene.Hou et al. (2011) reported a novel variant (T2643C and 2677AGTTCCCC deletion) in the intron 2 KISS1 gene in Chinese goat breeds.The quantified output indicated that TC genotype was correlated with high LS in Chinese goat breeds.Maitra et al. (2014) discovered a SNP at g.2803A>G in Indian goats.Contrary to the previous study, the present research could not find a mutation and an indel mutation at g.2643T>C, g.2677Del_AGTTCCCC, and g.2803A>G in Indonesian native goat breeds.This proposes a specific mutation characteristic of KISS1 gene within goat breeds in dissimilar environments related to the goat adaptation.
The mutation and deletion in the intron 1 KISS1 gene could influence the stability of mRNA both directly or indirectly.As a result, they impact protein production (Hou et al., 2011).SNP4 and SNP6 had no valuable effect on the prolific trait.Our result was consistent with the prior study in Chinese goat breeds, presenting that  these loci may not influence prolific traits.Furthermore, SNP11 and SNP12 locus significantly predisposed LS (p<0.05).Similar to this finding, SNP11 and SNP12 significantly influenced LS (p<0.05) in Chinese goat breeds.Three genotypes (TT, TC, and CC) have been identified at SNP11 in Indonesian goats.Contrary to this finding, only two genotypes are found (TT and TC) in Chinese goats.Nevertheless, the superior TC that affected LS was only found in Saanen goats, similar to Indonesian goats (An et al., 2013).The various distributions of genotype diversity are related to the different selection objectives, which resulted in the different breeds (Martínez-Royo et al., 2012).
In the recent research, we identified four novel SNPs of KISS1 gene intron 1 and two of them are in prominent association with reproductive traits in Indonesian native goat breeds.AA genotype at SNP9, GA genotype at SNP10 significantly affected LS in all Indonesian native goat breeds (p<0.01).Meanwhile, the genotypes of SNP7, SNP8, and SNP12 did not have the same correlation with LS in different goat breeds.Interestingly, AG genotype at SNP7, CC genotype at SNP8, and AA genotype at SNP12 impressed higher LS only in KC and KJ breeds.Additionally, CC genotype at SNP8 at SD goats has governed the same LS average as KC (2.67±0.23)and is higher than KJ (2.33±0.25).SNP13 was only found in KJ goats, and CC genotype affected the LS significantly (p<0.05).Conversely, TT genotype in Chinese goat breeds had greater LS than those CC in the fourth parity (An et al., 2013).
Otherwise, in SD goats, GG seems superior to AG at SNP7 and monomorphic at SNP8 and SNP12.However, the average LS (2.67±0.23) is likely as high as two other breeds.TT genotype at SNP14 had a larger LS than TA.The population is not under HWE at g.1960-1977DelTATGAGTAGCCCCCCTGC, SNP8, SNP9, SNP11, and SNP12.Furthermore, TT genotype was the most common genotype at 2601 st position (0.75).Selection in the screen population for higher LS and milk production has resulted in the present situation.Moreover, these mutations may be in linkage disequilibrium with other causative SNPs both in the same gene or other neighboring gene because the polymorphism is not the direct cause for the phenotypic variation.
Comparing genotype frequencies between three goat breeds shows that SD had obvious differences from KC and KJ.Following Hou et al. (2011), Boer goats had marked differently from Saanen and Guanzhong goats based on genotype frequencies in the intron 1 KISS1 gene.Genotypes could define the polymorphisms that were absurdly correlated with goat breeds by statistical discrepancies (χ2 test).
Goats with ID genotype at g.1960-1977DelTAT-GAGTAGCCCCCCTGC had significantly higher LS than II and DD genotype.Wild-type homozygous genotypes had higher LS (p<0.05) in most SNPs except at SNP10 and SNP11, both in the first and the second parities.This denotes homozygous mutant variants have a negative effect on prolific traits.According to this finding, heterozygous mutations could escalate ovulation rate and LS, whereas, in homozygous circumstances, the effects may be contradictory, from sterility to hyperprolificacy (Majd et al., 2019).
The indel diversity was lesser than SNPs (Table 2).This locus can be identified in three Indonesian native goat breeds and exhibited a high varied D allelic frequency (56%-94%).The ID genotype can be judged as predominant genotypes in Indonesian native goat breeds.Similar to this finding, ID genotype of SPEF2 gene and CDC25A gene was encountered to correlate with higher LS than other genotypes (Chen et al., 2019;Yang et al., 2021).Contrary to this result, homozygous genotypes significantly had the highest LS are DD genotype of KISS1 gene in Chinese goat breeds, II genotype of casein alpha S1 gene, and II genotypes of KDM6A gene (Cao et al., 2010;Wang et al., 2018;Cui et al., 2018).The important effect of indel position might be related to the efficiency of transcription factors binding there from altering the certain gene function and affecting firstborn LS (Hui et al., 2020).
The mutations are more in intronic regions than exon regions, confirming the number of SNPs (Evans et al., 2014).It has been recommended that introns might withstand mutational disruption, probably by being mutation buffer, cover the exons and regulators of gene expression, resources of protein diversity, and controlling phenotypic variations (Jo et al., 2015;Wang et al., 2018).An intron is a noncoding sequence that has enhanced the function through evolution to maintain its maintenance (Parenteau et al., 2019).The importance of the SNPs identified in introns regions cannot be rejected as it has been found that mutations may change gene expression by impressing splicing accuracy/efficiency or codon biases (Cartegni et al., 2002).
Intron mutations may influence the mRNA level of genes by turning the methylation of CpG sites and changing binding transcription factors (Van Laere et al., 2003).Intron could result in existing microRNA (miRNA) and non-coding RNA binding sites (Hui et al., 2020) and may alter spatial configuration and position of the chromosome, thereby influencing the binding of augmenter to transcription factors and promoters (Spielmann et al., 2018).Sequences in the non-coding region can affect the mechanism of mRNA deadenylation and degradation (Xu et al., 1997;Nackley et al., 2006) and consequently influence the amount of protein produced.Thus, intronic mutations should be considered for association with the functional traits (Maitra et al., 2014).
Point mutations at transcription factor binding sites cause changes in transcriptional regulation, resulting in phenotypic variation between breeds and species (Wray et al., 2007).According to statistical analysis (Table 7), the A allele at SNP9 strongly correlates with litter size.Therefore, the polymorphism at SNP9 may affect the binding affinity of TF.The mRNA profiling analysis on wild type animals revealed that NFIC is a more potent gene activator than a repressor, suggesting that NFIC may directly activate gene expression (Pjanic et al., 2011), which indicated the wild type A allele of SNP9 could be a significant genetic regulator of KISS1 gene expression through modulating kisspeptin that can affect GnRH release, and as a result, increasing larger lit-ter size.In addition, the g.1384G allele of the KISS1 gene may disrupt a GZF1 transcription factor binding site, thus reducing litter size (An et al., 2015 b ), and g.2540T had putative binding sites for the androgen receptor, which plays an important role in the signaling pathway involved in ovulation rate increase (Jeet et al., 2022).
HWE is utilized to appraise current populations' equilibrium and genotype frequencies and forecast allele and genotype frequencies in the next generations and earlier generations (Douhovnikoff et al., 2016).The χ2 values (Table 3) showed that most populations were not under HWE at five loci (g.1960-1977DelTATGAG-TAGCCCCCCTGC, SNP8, SNP9, SNP11, andSNP12).KC goats were in agreement with HWE at all loci.KJ goats were not under HWE at SNP8, SNP9, and SNP12, and neither were SD goats at g.1960-1977DelTATGAG-TAGCCCCCCTGC, SNP8, SNP9, SNP11, and SNP12 (χ2>3.841).Surprisingly, the disequilibrium sites in SD goats are the same as disequilibrium at the overall sites.A population in the covenant of HWE shows the lack of physical and natural selection, genetic drift, and equal genotypic distribution between generations (El-Komy et al., 2021).Moreover, major artificial selection of goat production traits might lead the population not to be under Hardy Weinberg Equilibrium (Yang et al., 2018;Dong & Du, 2011).The output of selection pressure caused higher dairy quantity and body weight than two other breeds, so SD goats are utilized as a dual-purpose goats.Parallel to this report, the Istrian Pramenka breed in Slovenia resulted from dairy capacity selection pressure (Starič et al., 2020).The selection program in KJ goats revealed a higher average body weight than KC goats (Kurnianto et al., 2013).Table 4 showed that the twin percentage on KJ goats is higher than KJ and SD goats; moreover, the single kid percentage is lower than two other breeds.This situation could be explained by more intensive economic traits selection in KJ and SD goats than KC goats.
Linkage disequilibrium (LD) is a strong method used to calculate correlation in many loci.The present study found a strong LD between SNPs associated with the prolific trait in Indonesian native goat breeds, which were SNP1 and SNP13, SNP8 and SNP9, SNP8 and SNP10, SNP8 and SNP12, SNP9 and SNP10, SNP9 and SNP12, SNP10 and SNP12 also SNP11 and SNP14, whereas SNP1 and SNP9 are novel variants of KISS1 gene.D' data ranged from 0.76 to 1, indicating very high linkage (Table 8).Loci that were linked firmly to profitable mutations can increase quickly in frequency by genetic hitchhiking (Waples et al., 2015) and interpreting their entity and coinheritance effect on the studied phenotypic traits (El-Komy et al., 2021).An intronic variation could be in perfect LD with known phenotypeassociated mutations (Nakaoka et al., 2016;Cui et al., 2018).LD leftover even after a few next generations of recombination and genome crosses, the advantages of these markers in large genetic selection of the pattern of segregating family are negligible (de Lima et al., 2020).
The present study analyzed the association of different genotypes and parities (from the first to the third parity).LS tended to incline in the third parities.Although, this trend was not found from the first parity to the second parity.The LS and puberty are governed by many factors such as nutritional, genetic, environmental, and maternal effects (de Lima et al., 2020).The effect of parity and year of kidding on LS was significant, but it varied depending on the loci combination.CC genotype at SNP13 is a marvelous genotype significantly correlated with larger LS in the first to the third parity (p<0.05).Unfortunately, it is found solely in KJ goat.Maitra et al. (2014) also found the disparity in LS and parity found in different genotypes.Starič et al. (2020) reported that genotype has different effects on different breeds, similar to the current study.These differences were correlated with genotype or allele frequencies and the correlation between reproductive performance and genotype diversity.This was also confirmed by Jeet et al. (2022) that the effect of the genotypes on reproductive functions could be influenced by the objective of the breed, environmental condition and/or the previous selection of the observed animals (Hernandez et al., 2005).SNP7, SNP11, and SNP12 demonstrated an incline of LS average from the first to the third parities.Although, the superior genotype correlated with higher LS and larger average LS per parity was not the same for those loci.Furthermore, these SNPs were in low LD (D'= 0.042, r= 0.00).Novel SNP1 and SNP14 significantly associated with LS (p<0.05),regrettably both SNPs were in low LD (D'= -1.00, r= 0.02).Therefore, these SNPs could not be used as a prominent genetic marker.
On the other hand, SNP8, SNP9, and SNP10 expressed that genotype at these loci influences the LS significantly at the first and the third parities, but not at the second parity.However, GA genotype at SNP10 showed an effect only in the first parity.Moreover, CC genotype at SNP8 and AA genotype at SNP9 showed an association with LS at the third parity, which was significantly higher than LS at the second parity and the first parity.Unfortunately, we were lack of samples at the third parity.Overall, CC genotype at SNP8, AA genotype at SNP9, and GA genotype at SNP10 were excellent genotypes associated with LS and parity significantly (p<0.05) in Indonesian native goat breeds.Additionally, these loci pairs (SNP8 and SNP9; SNP8 and SNP10; SNP9 and SNP10) were in remarkably strong linkage.(D'= 1 respectively, r= 1;0.58; 0.58 each, χ2= 23.000; 13.382; 13.382).In contrast, the interaction effect between parity and genotype was insignificant on LS in Black Bengal goats (Maitra et al., 2014).Selection for LS within the breed is suggested as a main genetic improvement strategy if only the environmental adaptation had been long-built and remains crucial, or if a 0.1-kid increase in LS is sufficient to increase production (Notter, 2012).
The reproductive traits are complex quantitative traits involving many genes, loci, and interactions.Therefore, it is necessary to analyze multiple genes or loci's combined effect on reproductive traits (An et al.,  2015 a ).

CONCLUSION
Fifteen variants have been recognized in the present research, of which four are novel SNPs (SNP1, SNP2, SNP3, and SNP9).SNP9 at intron 1 of KISS1 gene was found to be in strong linkage with SNP8 and SNP10, and identified to have a significant association with LS and parity.Furthermore, the does with H2 haplotype (CCATAGCGCAACGT) had higher LS than those other haplotypes.The NFIC at wild type A allele of SNP9 might act as an activator of gene expression.

Table 1 .
SNPs detected in intron 1 of KISS1 gene in Indonesian native goat breeds Note: SNP position was based on ENSCHIT00000037363.1 (known variants in this transcript) or NCBI (GenBank Accession Number GU142847.1)for the unknown variants in this transcript).

Table 2 .
Allele distribution and genotype frequencies of the SNPs analyzed in the KISS1 intron 1 of Indonesian goats Note: KC= Kacang goats; KJ= Kejobong goats; SD= Senduro goats; Overall= allele and genotype frequencies of three goat breeds; n= number of samples.The genotyped individuals are represented by the numbers in the brackets.toHWE,aswas indicated by χ2 calculation was greater than χ2 in the Table3.The expected heterozygosity (He) ranged from 0.04 to 0.45.KC goats were in HWE in all SNPs sites.KJ goats were not under HWE at SNP8, SNP9, and SNP12 (χ 2 >3.841), while SD goats were not appropriate in HWE at g.19601977DelTATGAGTAGCC CCCCTGC, SNP8, SNP9, SNP11, and SNP12 (χ 2 >3.841).

Table 3 .
Expected heterozygosity and Chi Square of SNPs determined at intron 1 of KISS1 gene : E= equilibrium in HWE, if χ2 value less than 3.841 (χ2 distribution table: 0.05; df= 1); D= disequilibrium in HWE if χ2 value more than 3.841; (-)= the χ2 value could not be calculated because the loci are monomorphic; Overall= the value of three Indonesian native goat breeds. Note

Table 4 .
The mean litter size in three Indonesian native goat breeds and frequency of birth type Note: Overall= mean of three goat breeds; SE= standard error.

Table 5 .
Means ± SE of litter size based on breeds, parities, and farms

Table 6 .
Least square means and standard error of litter size of parities 1 to 3 in Indonesian native goat breeds Note: Means with different superscripts at the same column and SNP differ significantly (p<0.05);N= number of samples; SE= standard error; NM= not measured (no samples available during investigation).

Table 7 .
Average litter size ± SE on different genotypes of KISS1 gene on Kacang, Kejobong, and Senduro goats Note: Means in the same column and same SNP followed by lower case letters vary significantly (Tukey test, 5%); N= number of sample for overall means; Overall= mean of LS of three goat breeds; The bold marked SNPs were included for further statistical analysis to ensure the association with prolific trait; SE= standard error; NM= Not measured (no sample available during investigation).

Table 9 .
Association of haplotype polymorphisms of KISS1 gene with litter size (p<0.0001)