Genomic Signatures of Selection Associated With Litter Size Trait in Jining Gray Goat

Abstract

Litter size (LS), an important economic trait in livestock, is so complicate that involves many aspects of reproduction, the underlying mechanism of which particularly in goat has always been scanty. To uncover the genetic basis of LS, the genomic sequence of Jining Gray goat groups (one famous breed for high prolificacy in China) with LS 1, 2, and 3 for firstborn was analyzed, obtaining 563.67 Gb sequence data and a total of 31,864,651 high-quality single nucleotide polymorphisms loci were identified. Particularly, the increased heterozygosity in higher LS groups, and large continuous homozygous segments associated with lower LS group had been uncovered. Through an integrated analysis of three popular methods for detecting selective sweeps (Fst, nucleotide diversity, and Tajima’s D statistic), 111 selected regions and 42 genes associated with LS were scanned genome wide. The candidate genes with highest selective signatures included KIT, KCNH7, and KMT2E in LS2 and PAK1, PRKAA1, and SMAD9 in LS3 group, respectively. Meanwhile, functional terms of programmed cell death involved in cell development and regulation of insulin receptor signaling pathway were mostly enriched with 42 candidate genes, which also included reproduction related terms of steroid metabolic process and cellular response to hormone stimulus. In conclusion, our study identified novel candidate genes involving in regulation of LS in goat, which expand our understanding of genetic fundament of reproductive ability, and the novel insights regarding to LS would be potentially applied to improve reproductive performance.

Introduction

Goat (Capra hircus), one of the most important livestock animals, has been domesticated from wild goat near Fertile Crescent in Middle East around 10,000 years ago (Naderi et al., 2008; Daly et al., 2018). Accompanied by the migration activity of human, the goat exhibited the strong adaptability to diverse environments, and gradually spread all over the world (Benjelloun et al., 2015). To meet human demand, many functional traits, such as dairy, meat, cashmere yield, and reproductive performance, all underwent a long period of natural selection and artificial intervention, ultimately, some specialized goat breeds with high performance of production were generated, which provided abundant genetic resources for modern domesticated goat.

Litter size (LS, numbers of lambing per ewe) is such a complex trait, which involves multiple biological process, including hormone secretion, follicular growth, ovulation, fertilization, embryo implantation, placenta and fetal development, and so on (Feng et al., 2015; Xu et al., 2018). Among which, the ovulation rate and LS are most important (Fogarty, 2009). Previously, more attentions have been attracted in reproductive trait study of sheep, four major prolificacy (that is LS) genes have been identified, that are bone morphogenetic protein receptor 1B (BMPR1B, FecB), bone morphogenetic protein 15 (BMP15, FecX), growth differentiation factor 9 (GDF9, FecG), and beta-1,4-N-acetyl-galactosaminyl transferase 2 (B4GALNT2, Fec^L) (Polley et al., 2010; Drouilhet et al., 2013; Våge et al., 2013; Abdoli et al., 2016). In particular, it has been proved an extra ovulation increased per estrus for one copy of BMP15 (0.4) and BMPR-1B (1.5) mutation (Davis, 2005; Abdoli et al., 2016; Talebi et al., 2018). Actually, the mutations of GDF9 and BMP15 genes had been verified to be strongly associated with LS trait in sheep (Demars et al., 2013; Våge et al., 2013). Recently, one research of genome-wide association studies (GWAS) in sheep had focused on LS trait, the high- and low-prolificacy sheep breeds were applied to explore the significant variants associated with LS, it revealed amount of candidate single nucleotide polymorphisms (SNPs), including BMPR1B, FBN1, and MMP2 in Wadi sheep, SMAD1 and CTNNB1 in Hu sheep, NCOA1 in Icelandic sheep, FLT1, NF1, PTGS2, and PLCB3 in Finn sheep, ESR2 in Romanov sheep, and ESR1, ETS1, FLI1, SPP1, and MMP15 in Texel sheep (Xu et al., 2018).

Currently, the genetic research on the characteristics of goat is lagging, and the limited information of LS trait in goat mostly derived from reports of sheep, that major in BMPR1B, BMP15, and GDF9 (An et al., 2015a, b; Ahlawat et al., 2016). Besides, there were few association analysis of single gene, for example, the SNPs of POU class 1 homeobox 1 (POU1F1) and KISS1 genes, insertion/deletion (indel) variants of lysine demethylase 6A (KDM6A) and β-Type platelet-derived growth factor receptor (PDGFRB), and transcript variants of catenin beta 1 (CTNNB1) had been reported to be associated with goat LS trait (An et al., 2013; Cui et al., 2018; Zhang X. et al., 2018; Yang et al., 2019; Zhu et al., 2019). Meanwhile, the genome-wide research on goat is also very lacking, only fewer reports showed some novel mutation sites affecting goat reproduction by whole-genome sequence (Guan et al., 2016; Lai et al., 2016). According to Guan et al. (2016), the reproduction-related genes regard to neurohypophysial hormone activity, including PAIP2B, CCDC64, and EPB41L5, had been identified. Also, Lai et al. (2016) revealed a series of candidate genes for fecundity (number of lambs) in Laoshan dairy goats, such as synonymous mutation: MAP3K12, DNMT3B, and SMAD2 and non-synonymous mutation: SETDB2 and CDH26. Moreover, PSEN2 gene, that locates at surrounding transcription start sites, was suggested to be related to goat fecundity, PRP1 and PRP6 genes were discovered to be linked to reproductive hormone with copy number variation (CNV) analysis (Zhang R.-Q. et al., 2018; Zhang et al., 2019).

Jining Gray goat, living in Shandong province of China, is one of the most outstanding local breeds in the world, because of its high fertility performance and it has shorter sexual maturity than other breeds (Miao et al., 2016; Su et al., 2018), thus the breed is the optimal animal model for detection of selection signatures associated with reproductive ability in goat. Several reports have investigated the reproductive performance focused on Jining Gray goat from different aspects (Chu et al., 2011; Shi et al., 2015; Miao et al., 2016; Su et al., 2018), but the variant research in genome-wide has not been reported. Actually, as one kind of low heritability traits in livestock, LS has been hard to improve with traditional selection. Presently, the technologies of genomic selection (GS) in livestock have become the most promising power to improve the traits with low heritability (Cm Dekkers, 2012; Jonas and Koning, 2015). Despite of under development, the implementation of GS in livestock has been effective for increasing the accuracy of estimated breeding values (EBVs), by using the genome-wide SNP in breeding programs, it shortens generation intervals and is more economical (VanRaden et al., 2009; Pryce and Daetwyler, 2012). Regard to this, the effective genetic variants associated with breeds advantage genotypes have to be uncovered, especially in goat.

To elucidate the genetic fundament of LS in goat, 40 Jining Gray goats, with LS of 1 (20 individuals), 2 (14 individuals), and 3 (6 individuals) for firstborn, were performed with the whole-genome sequence at ∼10× coverage in this study. In the result, amounts of whole-genome variants, including SNPs and indel, had been discovered. Meanwhile, the selective sweeps of different LS groups were contributed to screen the regions with strong signatures of selection, and the most candidate genes had been identified within Jining Gray goat population, that provided abundant candidate loci for reproduction improvement in domestic goat.

Materials and Methods

Ethics Statement

All goats in this study were keep feeding in Caoxian Zhengdao Jining Gray goat’ farm in Heze, Shandong, China. The experiment procedures were approved by the Ethics Committee of Northwest A&F University.

Animals

The selected goats were 2∼3 years old female individuals, and the natural mating was arranged with normal reproductive male goats. And the LS information was derived from record of Caoxian Zhengdao Jining Gray goat farm. Total 40 Jining Gray goats were divided into three groups with LS of 1 (20 individuals), 2 (14 individuals), and 3 (6 individuals) for firstborn.

Samples Collection, DNA Extraction, and Whole-Genome Sequence

Ear tissues of selected goats were collected in farm, transferred in sterile centrifuge tube, and stored in library at −80°C. For sequence, the samples were transferred with dry ice to the Beijing Novogene Bioinformatics Technology Co., Ltd. The DNA sample of each individual was extracted and measured following standard experimental procedures to guarantee their quantity. Subsequently, the DNA samples were interrupted randomly to produce fragments of 350 bp, and the library was constructed with a series of steps of end repair, ployA tailed, purification, PCR amplification, and so on. The whole-genomic sequence with paired-end 150 bp strategy was performed by the Illumina NovaSeq 6000 platform (Illumina Inc., San Diego, CA, United States).

Genomic Data Processing and Variants Calling

With the sequence data from company, the Trimmomatic software (version 0.38) (Bolger et al., 2014) was applied to remove adaptors and low-quality reads, and FastQC reports (version 0.11.6) (Andrews, 2010) were generated to check data quality (for statistics information, see Supplementary Table S2). Through BWA software (version 0.7.17) (Li and Durbin, 2009), the clean data were aligned to goat reference genome¹ (ARS1) (Bickhart et al., 2017), and the bam and sorted bam files were generated with samtools (version 1.9) (Li et al., 2009). Then, the Genome Analysis Toolkit (GATK, version 4.1.3.0) (McKenna et al., 2010) pipeline was executed with diverse modules for variants calling: briefly, it was prepared with duplicates fix of sorted bam files, then “HaplotypeCaller” module produced the raw variants of gVCF files, and “CombineGVCFs” was carried to merge all raw results as one gVCF file, then “GenotypeGVCFs” module had been operated to consolidate the calling variants to generate VCF file (including SNPs and indels). Finally, the preliminary SNPs and indels were produced with the VCF file by “SelectVariants” mode.

Quality Control and Annotation for Detected SNPs

To guarantee the quality, the raw SNPs were filtered by a series of thresholds with “VariantFiltration” mode of GATK. That is: “QD (Quality score by depth) >2.0 & MQ (Mapping quality score) >40.0 & FS (Phred-scaled p-value using Fisher’s exact test) <60.0 & ReadPosRankSum (Rank Sum Test for relative positioning of reference and alternative alleles within reads) > −8.0 & MQRankSum (Rank sum test for mapping quality of reference and alternative within reads) > −12.5,” which had been widely used to generate the high-quality SNPs for subsequent analysis. To identify the variants function, ANNOVAR software was utilized to annotate the filtered SNPs (version 2018-04-16) (Wang et al., 2010). The SNPs were displayed with circos plot with shinyCircos² (Yu Y. et al., 2017) and density plot by CMplot package³ in R environment (v 3.6.0).

Population Genetics and LD Analysis

For principal component analysis (PCA), the software of plink (v1.90b4) (Purcell et al., 2007) was performed to prepare the input file (.bed, .bim, .fam) and conduct simple filtering (–geno 0.05), with autosomal SNPs of 95% genotyping rate, gcta64 (v 1.26.0) (Yang et al., 2011) was applied to calculate genetic relationship matrix (–make-grm) and principal component (–pca). And plot was generated with python in local. The neighbor-joining (NJ) phylogenetic tree was constructed by MEGA (v 10.0.5) (Saitou and Nei, 1987; Kumar et al., 2018) software based on the genetic distance-matrix calculating by plink and plotted with iTOL (v 4.4.2) (Letunic and Bork, 2019). Also, the linkage disequilibrium (LD) decay of r² (a measure of LD, that is squared coefficient of correlation (r) between loci) was carried out with PopLDdecay with general parameter (Zhang C. et al., 2018).

Runs of Homozygosity

The runs of homozygosity (ROH) value was estimated with autosomal SNPs of each individual by plink (v1.90b4). In our research, six criteria were set to define a ROH: a window size of 50 SNPs for scanning; the minimum length of ROH was 200 kb; the required minimum density was 50 kb/SNP; the maximum gap between consecutive homozygous SNPs was 500 kb length; a scanning ROH can contain at most one heterozygous and five missing genotypes; threshold of all scanning windows containing SNP was 0.05. As for assessment, ROH were divided into five categories based on length: 0–0.5 Mb, 0.5–1 Mb, 1–2 Mb, 2–4 Mb, and >4 Mb. All the parameter setting and classification for assessment had referenced other published reports (Peripolli et al., 2017; Onzima et al., 2018; Xu et al., 2019).

Detection of Selective Signatures

In preliminary, the values of pairwise Fst, heterozygosity, allele frequency, nucleotide diversity and Tajima’s D were calculated by vcftools with parameter of “--weir-fst-pop,” “--het,” “--freq,” “--window-pi,” and “--TajimaD” (v 0.1.16) (Danecek et al., 2011). For heterozygosity rate, it was calculated by the formula of (N_SITES-O(HOM))/N_SITES, O(HOM) represents the observed number of homozygous sites of each individual, N_SITES is the total number sites of each individual. With comparisons between diverse LS groups, the selective sweeps analysis was performed according the genetic differentiation (Fst), nucleotide diversity, and Tajima’s D values, which were computed with 100 kb window size and 50 kb step size. And the per window Fst was standardized to Z-transformed Fst [Z(Fst)] based on the formula of “per window Z(Fst) = (per window Fst − mean Fst of all windows)/standard deviation of all windows Fst.” The threshold of Z(Fst) value = 5 was applied to detect the selective regions. Moreover, the highly selected regions were supported by Tajima’s D values of different LS groups and pairwise nucleotide diversity (π). Manhattan plots of Z(Fst) value for comparisons between groups were made with the R package: qqman (Turner, 2014).

Functional Annotation for Candidate Genes

To uncover the biological functions of candidate genes, we performed the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis with Metascape (Zhou et al., 2019), a web-based portal for functional enrichment.⁴ Due to the inadequate database of enrichment analysis in goat, the goat gene symbol IDs were converted into mouse homologous gene symbol IDs through Ensembl website.

Results

Sequence and Variants Detection of Jining Gray Goat

In the present study, we sequenced 40 Jining Gray goats, which were divided into three groups with LS values of 1, 2, and 3 for firstborn. The total 563.67 Gb raw data was generated with average sequence depth >8 of each individual sample (Supplementary Tables S1, S2), and mean depth of LS groups mainly ranged from 9 to 11 (Supplementary Figure S1A). With high mapping rate (over 99%) to reference genome, we called 32,072,110 SNPs, after filtered, we got 31,864,651 SNPs with high quality (Table 1). Among which, the substitution mutations with largest number were C > T and G > A, secondly, A > G and T > C followed (Table 1). Moreover, the densities of SNPs and indels spreading on chromosomes were displayed with circos plot (Figure 1A), which indicated the maximum density of SNPs and indels occurred on chromosomes 18 and 23, respectively.

FIGURE 1

Moreover, the high-quality SNPs in autosome reached 28,926,281 Supplementary Figure S1B), and the SNPs number within 1 Mb window size was mainly around 8,000 and high-density SNPs (more than 16,000 in each window) largely distributed at the ends of chromosomes (Supplementary Figure S1C). With annotation, we got 9,447,885 SNPs variant results (29.65% of total) and 8,878,492 at autosome (Supplementary Figure S1B), which meant about two-thirds were novel SNPs in goat. These SNPs were partitioned into intergenic (6,343,642; 67.14%), intronic (2,480,758; 26.26%), exonic (64,074; 0.68%) and other gene regulatory regions (Figure 1B and Table 1). For SNPs in exonic regions, it showed the overwhelming majorities were synonymous and non-synonymous, and others including stopgain, stoploss, and unknown were <1% (Figure 1C and Table 1).

Divergent Genetic Differentiation of Jining Gray Goat LS Groups

To disclose the genetic differences between groups with different LS values, the SNPs of diverse LS groups were investigated. As showed in Figure 2A, there were more SNPs in LS2 than LS1 group, and LS3 group was medium. Similarly, the heterozygous/homozygous (het/hom) ratios in LS2 and LS3 were a little bit higher than that of LS1 group (Figure 2B), as well as the heterozygosity rates of different groups (Figure 2C). ROH analysis further revealed large ROHs of >4 Mb were most associated with lower LS group (Figure 2D), that supported the increase of heterozygosity rates in high LS groups. Moreover, we performed PCA based on SNPs in autosomes (Figure 2E), indicating these samples had no obvious structure regarding the LS groups, which was consistent with the NJ tree result (Supplementary Figure S2A). Further, the LD analysis with LD decay of r² showed there was no apparent difference among the three LS groups (Supplementary Figure S2B), which was likely because these individuals come from the same large population.

FIGURE 2

Signatures of Selection With Sweeps Analysis

To detect signatures of selection associated with LS, the genetic differentiation was estimated with the fixation index (Fst) in groups of Jining Gray goat (Supplementary Table S3). Firstly, genome-wide SNPs in autosome were utilized to calculate Fst value of LS2 vs LS1 groups, which was performed with 100 kb window size and 50 kb step size. As previous reports, the regions with high Z-transformed Fst values [Z(Fst) > 5] were concentrated (Axelsson et al., 2013; Gou et al., 2014). Here, total 67 windows [Z(Fst) > 5] had been strongly selected (Figure 3A), where KCNH7, KMT2E, and KIT genes located, other genes, consist of CNOT9, SULT4A1, and SERPINB10, underwent subordinate selection (Table 2 and Supplementary Table S4). Especially, there were some novel genes in goat, part of which were generated with the homologous in mouse, i.e., Akr1c21 and Serpinb2 (Figure 3A). Genes, KIT, KCNH7, and CNOT9, were verified with Tajima’s D values and π (pairwise LS2 vs LS1) (Figure 3B). Allele frequencies for these loci are shown in Supplementary Figure S2C. Moreover, the non-coding region of KIT and KCNH7 also possessed many mutations (Supplementary Figures S3A,B). Among, KIT gene was most likely the candidate gene in LS2 group of Jining Gray goat.

FIGURE 3

Further, the selective sweep analysis also was executed in LS3 vs LS1 groups (Figure 4A and Table 2). Here, we obtained 44 regions with highest Z(Fst) value above the genome-wide line (Supplementary Table S3), which located at chromosome 9, 20, and 29, and the selection genes of SMAD9, CARD6 and PRKAA1, and PAK1 were annotated, respectively. Also, the genes at a sub-optimal level included PARD3B, SORL1, and PSMA6 (Table 2 and Supplementary Table S4), especially. Several novel LS genes in goat were mapped in mouse orthologs genes, such as Pctp and Gdpd4. Moreover, these genes in selective regions of LS3 group, PAK1 and Gdpd4, CARD6 and PRKAA1, and PSMA6, were partly confirmed by Tajima’s D and π values compare to LS1 group (Figure 4B), that also had been reflected by the changes of allele frequency at mutation loci (Supplementary Figure S2D). In addition, we discovered that PAK1 gene had several mutations in non-coding region (Supplementary Figure S3C). These results indicted PAK1 was mostly candidate gene in LS3 group of Jining Gray goat.

FIGURE 4

Functional Enrichment of Candidate Genes

To disclose the function of candidate genes associated with LS trait, the GO enrichment analysis of annotated genes in regions with Z(Fst) > 5 both LS2 vs LS1 and LS3 vs LS1 groups was performed (Table 3). The most top of non-redundant enrichment GO terms were “programmed cell death involved in cell development” and “regulation of insulin receptor signaling pathway” (Figure 5A). Of note, enrichment map of all GO lists showed “steroid metabolic process,” “cellular response to hormone stimulus,” and “developmental process involved in reproduction” were also enriched with more gene hits (marked with numbers) (Figure 5B and Supplementary Table S5), all these processes were likely evidenced the genotype of the high LS of firstborn in Jining Gray goat. Meanwhile, the enriched pathway of “tight junction” might also serve as other important roles in high LS (Figure 5C and Table 4).

FIGURE 5

Discussion

In this study, we provided the comprehensive evidence of genetic background of high fecundity in goat. Briefly, the whole-genomic sequence was performed with groups of Jining Gray goat (one famous Chinese local breed with mean LS 2.94) with LS 1, 2, and 3 for first-born and we got a plenty of genome-wide variants. Most of genomic SNPs in Jining Gray goat were identified for the first time in our study. Importantly, signatures of selection among LS groups were uncovered with selective sweeps of Fst, π, and Tajima’s D. We identified 42 candidate genes associated with LS trait, which was likely contribute for GS in goat to improve the reproductive performance.

As previously reported, follicular maturation and ovulation rate are the most important factors affecting LS, increasing the possibility of higher fertilized rate and bigger litters (Fogarty, 2009). It has been suggested that the members of transforming growth factor-β (TGF-β) superfamily, including BMPR-1B, BMP15, and GDF9, affected the ovulation and LS (Demars et al., 2013). Also, numerous studies had explored the polymorphisms of these genes in sheep and goat breeds (Ahlawat et al., 2016; Talebi et al., 2018). Here, we identified large amounts of SNPs in Jining Gray goat, and more than half were novel sites, all the abundant variants extended genetic resource in goat for fecundity traits associated research. Of note, there were no distinct grouped patterns in PCA and NJ tree results, as well as LD decay analysis, which might cause by the similar genetic background within population. Recently, ROH in animals present adamant evidence of individual relatedness, inbreeding, and also selection pressure (Peripolli et al., 2017; Ceballos et al., 2018). Our results showed that large continuous homozygous segments were mostly associated with lower LS group, which displayed consistent with increase of heterozygosity rate in higher LS groups, these results indicated the difference of heterozygosity in diverse LS groups was likely related to fitness and inbreeding.

For genomic research, it has been established that GWAS is the most powerful tools for detecting the genetics basis of complex traits (McCarthy et al., 2008; Stranger et al., 2011). Accordingly, an adequate statistical power of GWAS require a much larger sample size, and the proper case-control design is relative necessary (Hong and Park, 2012). However, these very strict requirements of the large sample size and detailed genealogy information prevent us to carry out the similar study in goat, regard to their availability. At the same time, the threshold character of LS trait had been reported as previously (Konstantinov et al., 1994; Pérez-Enciso et al., 1995; Janssens et al., 2004), and the extreme phenotypic study design is also much more powerful to detect the variants associated trait selection (Peloso et al., 2016), thus, the sweeps analysis of pairwise value between LS groups were conducted.

In our study, a total of 42 candidate genes associated with LS trait were identified by selective sweeps of three groups within Jining Gray goat. In LS2 group, the genes of KCNH7, KMT2E, and KIT possessed most strong selective signatures. KIT, is known to be involved in oocyte growth and follicular development in mammals, deficient of KIT would lead amounts of oocytes lost (Moniruzzaman et al., 2007), and evidence showed mutations of KIT or KITLG blocked the interaction of oocyte and granulosa cell, and caused female infertile (Horie et al., 1991; Driancourt et al., 2000; Thomas et al., 2008). Moreover, KITLG has been previously identified to be associated with LS trait in goat (An et al., 2015a, b), such evidence suggested KIT may be the most promising gene for large LS. Another, KCNH7 also had been found to be significantly associated with animal reproductive traits (Fan et al., 2017). Besides, KMT2E, termed as MLL5, had been demonstrated to be required for normal spermatogenesis, its knockout mouse showed a post meiotic phenotype and infertile (Heuser et al., 2009; Yap et al., 2011).

Besides, in LS3 group, the genes of SMAD9, CARD6, PRKAA1, and PAK1 possessed most strong selective signatures. Among which, SMAD9 (also known as SMAD8), a member of the SMAD family, involved signal transduction of TGF-β superfamily and inhibited BMP activity (Chang et al., 2002; Tsukamoto et al., 2014), it had already been suggested to control follicular selection via balancing luteinizing hormone receptor (LHR) transcription (Yu D. et al., 2017). Similarly, SMAD2 in goat and SMAD1 in sheep also were detected to be associated with LS trait (Lai et al., 2016; Xu et al., 2018). In addition, PAK1, a serine/threonine kinase, had been reported to mediate estradiol-negative feedback actions in the reproductive axis, which seems to play important roles in non-classical estrogen receptor α (ERα) signaling (Zhao et al., 2009). PRKAA1, protein kinase AMP-activated catalytic subunit alpha 1, served as an energy sensor that maintains energy homeostasis (Kahn et al., 2005; Krishan et al., 2014). And another gene, PARD3B with lower selection intensity here, had been identified with goat productive traits in another report (Guan et al., 2016).

Especially, our results suggested regions detected with Fst were not supported by Tajima’s D or π. These contradictory results may be caused by differences in the underlying theory among methods. Firstly, Fst, is developed for population differentiation, it compares the variance of allele frequencies within and between populations, and is suggestive of directional selection (Holsinger and Weir, 2009). Another, Tajima’s D is calculated based on frequency spectrum, when selection occurs, a population wide reduction in the genetic diversity follows, and many surrounding sites near the selected variant segregate at low frequencies, then the frequency spectrum shifts back to baseline over thousands of generations (Tajima, 1989). Conversely, population differentiation-based approaches (that is Fst) can detect more types of selection, including classic sweeps, sweeps on standing variants, and negative selection (Vitti et al., 2013). Thus, we believed the result of Fst was potentially acceptable, although they might be not suggested by other values at secondary level of significance. Moreover, the inconsistencies between results of LS2 vs LS1 and LS3 vs LS1 attracted more attentions. Here we may hypothesize the LS trait is not only quantitative, but also a threshold trait, it means different major genes and more minor genes when LS increase. In our result, there were no shared candidate genes between in LS2 and LS3 groups, it seemed to indicate the character of threshold trait of LS, in particular, the allele frequency of PAK1 steadily increased from LS 1 to 3, which suggested a quantitative trait of LS, although it only showed strong selective signal in LS3 group. However, the underlying mechanism needs further investigation.

Further, the functional enrichment analysis suggested the processes of programmed cell death involved in cell development, cellular response to hormone stimulus, steroid metabolic process and regulation of insulin receptor signaling pathway might play important roles in the regulation of LS trait. Notably, the steroid production is critical for body development and function through steroid hormones, consist of androgen, estradiol, progesterone, glucocorticoid, and aldosterone (Jamnongjit and Hammes, 2006; Ruiz-Cortés, 2012). There were six genes enriched in steroid metabolic process (Supplementary Table S5), KIT and PRKAA1 both with high selective signals that exhibited reproduction or energy-related regulation. In another aspect, programmed cell death plays a fundamental role in animal development, and plays a prominent role in development of fetal ovaries and postnatal ovarian cycle (Pru and Tilly, 2001; Fuchs and Steller, 2011). Moreover, the insulin signal that interacted with steroid metabolic process and cellular response to hormone stimulus was believed to be closely connected with folliculogenesis and ovarian function (Kezele et al., 2002; Dupont and Scaramuzzi, 2016).

Collectively, this research identified numerous novel variants in Jining Gray goat. Some candidate genes showed strongest differentiation signals, such as KIT, PAK1, PRKKA1, SMAD9, KCNH7, and KMT2E. Especially, KIT and PAK1 genes are mostly suggested to be applied to improve the reproductive performance of goats. Meanwhile, the steroid metabolic process and cellular hormone response enriched most genes at selective regions, which would help to explain the difference of LS trait in Jining Gray goat. With regard to application, we believe two steps should be considered: (1) GWAS with large sample size in goat to verify the accuracy of candidate genes; (2) establishment of the core breeding population based on genetic marker genes and reference population expansion. Accordingly, GS is a promising development in livestock breeding programs (Jonas and Koning, 2015), whose implementation of is based on the genome-wide markers of interested trait to conduct genomic EBVs, with daughter assessment, the selected individuals are arranged to reference population. Above of all, our findings will increase the genetic cognition of goat LS trait, contribute to understand the changes of genomic phenotypes under the adaption and artificial breeding process, which would benefit the improvement of reproduction performance in goat.

References

• 1

  AbdoliR.ZamaniP.MirhoseiniS.Ghavi Hossein-ZadehN.NadriS. (2016). A review on prolificacy genes in sheep.Reprod. Domest. Anim.51631–637. 10.1111/rda.12733

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 2

  AhlawatS.SharmaR.RoyM.MandakmaleS.PrakashV.TantiaM. (2016). Genotyping of novel SNPs in BMPR1B, BMP15, and GDF9 genes for association with prolificacy in seven Indian goat breeds.Anim. Biotechnol.27199–207. 10.1080/10495398.2016.1167706

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 3

  AnX.HouJ.GaoT.LeiY.SongY.WangJ.et al (2015a). Association analysis between variants in KITLG gene and litter size in goats.Gene558126–130. 10.1016/j.gene.2014.12.058

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 4

  AnX.HouJ.LeiY.GaoT.SongY.WangJ.et al (2015b). Two mutations in the 5-flanking region of the KITLG gene are associated with litter size of dairy goats.Anim. Genet.46308–311. 10.1111/age.12277

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 5

  AnX.MaT.HouJ.FangF.HanP.YanY.et al (2013). Association analysis between variants in KISS1 gene and litter size in goats.BMC Genet.14:63. 10.1186/1471-2156-14-63

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 6

  AndrewsS. (2010). FastQC: a Quality Control Tool for High Throughput Sequence Data.Cambridge: Babraham Institute.

  • Google Scholar

• 7

  AxelssonE.RatnakumarA.ArendtM.-L.MaqboolK.WebsterM. T.PerloskiM.et al (2013). The genomic signature of dog domestication reveals adaptation to a starch-rich diet.Nature495360–364. 10.1038/nature11837

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 8

  BenjellounB.AlbertoF. J.StreeterI.BoyerF.CoissacE.StuckiS.et al (2015). Characterizing neutral genomic diversity and selection signatures in indigenous populations of Moroccan goats (Capra hircus) using WGS data.Front. Genet.6:107. 10.3389/fgene.2015.00107

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 9

  BickhartD. M.RosenB. D.KorenS.SayreB. L.HastieA. R.ChanS.et al (2017). Single-molecule sequencing and chromatin conformation capture enable de novo reference assembly of the domestic goat genome.Nature Genet.49643–650. 10.1038/ng.3802

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 10

  BolgerA. M.LohseM.UsadelB. (2014). Trimmomatic: a flexible trimmer for Illumina sequence data.Bioinformatics302114–2120. 10.1093/bioinformatics/btu170

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 11

  CeballosF. C.JoshiP. K.ClarkD. W.RamsayM.WilsonJ. F. (2018). Runs of homozygosity: windows into population history and trait architecture.Nat. Rev. Genet.19220–234. 10.1038/nrg.2017.109

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 12

  ChangH.BrownC. W.MatzukM. M. (2002). Genetic analysis of the mammalian transforming growth factor-β superfamily.Endocr. Rev.23787–823. 10.1210/er.2002-0003

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 13

  ChuM.LuL.FengT.DiR.CaoG.WangP.et al (2011). Polymorphism of bone morphogenetic protein 4 gene and its relationship with litter size of Jining Grey goats.Mol. Biol. Rep.384315–4320. 10.1007/s11033-010-0556-6

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 14

  Cm DekkersJ. (2012). Application of genomics tools to animal breeding.Curr. Genomics13207–212. 10.2174/138920212800543057

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 15

  CuiY.YanH.WangK.XuH.ZhangX.ZhuH.et al (2018). Insertion/deletion within the KDM6A gene is significantly associated with litter size in goat.Front. Genet.9:91. 10.3389/fgene.2018.00091

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 16

  DalyK. G.DelserP. M.MullinV. E.ScheuA.MattiangeliV.TeasdaleM. D.et al (2018). Ancient goat genomes reveal mosaic domestication in the Fertile Crescent.Science36185–88. 10.1126/science.aas9411

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 17

  DanecekP.AutonA.AbecasisG.AlbersC. A.BanksE.DepristoM. A.et al (2011). The variant call format and VCFtools.Bioinformatics272156–2158. 10.1093/bioinformatics/btr330

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 18

  DavisG. H. (2005). Major genes affecting ovulation rate in sheep.Genet. Sel. Evol.37 (Suppl. 1), S11–S23.

  • Google Scholar

• 19

  DemarsJ.FabreS.SarryJ.RossettiR.GilbertH.PersaniL.et al (2013). Genome-wide association studies identify two novel BMP15 mutations responsible for an atypical hyperprolificacy phenotype in sheep.PLoS Genet.9:e1003482. 10.1371/journal.pgen.1003482

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 20

  DriancourtM.-A.ReynaudK.CortvrindtR.SmitzJ. (2000). Roles of KIT and KIT LIGAND in ovarian function.Rev. Reprod.5143–152. 10.1530/ror.0.0050143

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 21

  DrouilhetL.MansanetC.SarryJ.TabetK.BardouP.WoloszynF.et al (2013). The highly prolific phenotype of Lacaune sheep is associated with an ectopic expression of the B4GALNT2 gene within the ovary.PLoS Genet.9:e1003809. 10.1371/journal.pgen.1003809

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 22

  DupontJ.ScaramuzziR. J. (2016). Insulin signalling and glucose transport in the ovary and ovarian function during the ovarian cycle.Biochem. J.4731483–1501. 10.1042/BCJ20160124

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 23

  FanQ.WuP.DaiG.ZhangG.ZhangT.XueQ.et al (2017). Identification of 19 loci for reproductive traits in a local Chinese chicken by genome-wide study.Genet. Mol. Res.16:gmr16019431. 10.4238/gmr16019431

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 24

  FengT.CaoG. L.ChuM. X.DiR.HuangD. W.LiuQ. Y.et al (2015). Identification and verification of differentially expressed genes in the caprine hypothalamic-pituitary-gonadal axis that are associated with litter size.Mol. Reprod. Dev.82132–138. 10.1002/mrd.22451

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 25

  FogartyN. (2009). A review of the effects of the Booroola gene (FecB) on sheep production.Small Rumin. Res.8575–84. 10.1016/j.smallrumres.2009.08.003

  • CrossRef
  • Google Scholar

• 26

  FuchsY.StellerH. (2011). Programmed cell death in animal development and disease.Cell147742–758. 10.1016/j.cell.2011.10.033

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 27

  GouX.WangZ.LiN.QiuF.XuZ.YanD.et al (2014). Whole-genome sequencing of six dog breeds from continuous altitudes reveals adaptation to high-altitude hypoxia.Genome Res.241308–1315. 10.1101/gr.171876.113

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 28

  GuanD.LuoN.TanX.ZhaoZ.HuangY.NaR.et al (2016). Scanning of selection signature provides a glimpse into important economic traits in goats (Capra hircus).Sci. Rep.6:36372. 10.1038/srep36372

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 29

  HeuserM.YapD. B.LeungM.De AlgaraT. R.TafechA.MckinneyS.et al (2009). Loss of MLL5 results in pleiotropic hematopoietic defects, reduced neutrophil immune function, and extreme sensitivity to DNA demethylation.Blood1131432–1443. 10.1182/blood-2008-06-162263

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 30

  HolsingerK. E.WeirB. S. (2009). Genetics in geographically structured populations: defining, estimating and interpreting F ST.Nat. Rev. Genet.10639–650. 10.1038/nrg2611

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 31

  HongE. P.ParkJ. W. (2012). Sample size and statistical power calculation in genetic association studies.Genomics Inform.10117–122.

  • Google Scholar

• 32

  HorieK.TakakuraK.TaiiS.NarimotoK.NodaY.NishikawaS.et al (1991). The expression of c-kit protein during oogenesis and early embryonic development.Biol. Reprod.45547–552. 10.1095/biolreprod45.4.547

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 33

  JamnongjitM.HammesS. R. (2006). Ovarian steroids: the good, the bad, and the signals that raise them.Cell Cycle51178–1183. 10.4161/cc.5.11.2803

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 34

  JanssensS.VandepitteW.BodinL. (2004). Genetic parameters for litter size in sheep: natural versus hormone-induced oestrus.Genet. Sel. Evol.36543–562.

  • Pubmed Abstract
  • Google Scholar

• 35

  JonasE.KoningD.-J. D. (2015). Genomic selection needs to be carefully assessed to meet specific requirements in livestock breeding programs.Front. Genet.6:49. 10.3389/fgene.2015.00049

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 36

  KahnB. B.AlquierT.CarlingD.HardieD. G. (2005). AMP-activated protein kinase: ancient energy gauge provides clues to modern understanding of metabolism.Cell Metab.115–25. 10.1016/j.cmet.2004.12.003

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 37

  KezeleP. R.NilssonE. E.SkinnerM. K. (2002). Insulin but not insulin-like growth factor-1 promotes the primordial to primary follicle transition.Mol. Cell. Endocrinol.19237–43. 10.1016/s0303-7207(02)00114-4

  • CrossRef
  • Google Scholar

• 38

  KonstantinovK.ErasmusG.Van WykJ. (1994). Evaluation of Dormer sires for litter size and lamb mortality using a threshold model.S. Afr. J. Anim. Sci.24119–121.

  • Google Scholar

• 39

  KrishanS.RichardsonD. R.SahniS. (2014). Amp kinase (prkaa1).J. Clin. Pathol.67758–763. 10.1136/jclinpath-2014-202422

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 40

  KumarS.StecherG.LiM.KnyazC.TamuraK. (2018). MEGA X: molecular evolutionary genetics analysis across computing platforms.Mol. Biol. Evol.351547–1549. 10.1093/molbev/msy096

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 41

  LaiF.-N.ZhaiH.-L.ChengM.MaJ.-Y.ChengS.-F.GeW.et al (2016). Whole-genome scanning for the litter size trait associated genes and SNPs under selection in dairy goat (Capra hircus).Sci. Rep.6:38096. 10.1038/srep38096

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 42

  LetunicI.BorkP. (2019). Interactive tree of life (iTOL) v4: recent updates and new developments.Nucleic Acids Res.47W256–W259. 10.1093/nar/gkz239

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 43

  LiH.DurbinR. (2009). Fast and accurate short read alignment with Burrows–Wheeler transform.Bioinformatics251754–1760. 10.1093/bioinformatics/btp324

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 44

  LiH.HandsakerB.WysokerA.FennellT.RuanJ.HomerN.et al (2009). The sequence alignment/map format and SAMtools.Bioinformatics252078–2079. 10.1093/bioinformatics/btp352

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 45

  McCarthyM. I.AbecasisG. R.CardonL. R.GoldsteinD. B.LittleJ.IoannidisJ. P.et al (2008). Genome-wide association studies for complex traits: consensus, uncertainty and challenges.Nat. Rev. Genet.9356–369. 10.1038/nrg2344

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 46

  McKennaA.HannaM.BanksE.SivachenkoA.CibulskisK.KernytskyA.et al (2010). The genome analysis toolkit: a MapReduce framework for analyzing next-generation DNA sequencing data.Genome Res.201297–1303. 10.1101/gr.107524.110

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 47

  MiaoX.LuoQ.ZhaoH.QinX. (2016). Genome-wide analysis of miRNAs in the ovaries of Jining Grey and Laiwu Black goats to explore the regulation of fecundity.Sci. Rep.6:37983. 10.1038/srep37983

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 48

  MoniruzzamanM.SakamakiK.AkazawaY.MiyanoT. (2007). Oocyte growth and follicular development in KIT-deficient Fas-knockout mice.Reproduction133117–125. 10.1530/rep-06-0161

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 49

  NaderiS.RezaeiH.-R.PompanonF.BlumM. G.NegriniR.NaghashH.-R.et al (2008). The goat domestication process inferred from large-scale mitochondrial DNA analysis of wild and domestic individuals.Proc. Natl. Acad. Sci. U.S.A.10517659–17664. 10.1073/pnas.0804782105

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 50

  OnzimaR. B.UpadhyayM. R.DoekesH. P.BritoL. F.BosseM.KanisE.et al (2018). Genome-wide characterization of selection signatures and runs of homozygosity in Ugandan goat breeds.Front. Genet.9:318. 10.3389/fgene.2018.00318

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 51

  PelosoG. M.RaderD. J.GabrielS.KathiresanS.DalyM. J.NealeB. M. (2016). Phenotypic extremes in rare variant study designs.Eur. J. Hum. Genet.24924–930. 10.1038/ejhg.2015.197

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 52

  Pérez-EncisoM.FoulleyJ. L.BodinL.ElsenJ. M.PoiveyJ. P. (1995). Genetic improvement of litter size in sheep. A comparison of selection methods.Genet. Sel. Evol.2743–61.

  • Google Scholar

• 53

  PeripolliE.MunariD.SilvaM.LimaA.IrgangR.BaldiF. (2017). Runs of homozygosity: current knowledge and applications in livestock.Anim. Genet.48255–271. 10.1111/age.12526

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 54

  PolleyS.DeS.BrahmaB.MukherjeeA.VineshP.BatabyalS.et al (2010). Polymorphism of BMPR1B, BMP15 and GDF9 fecundity genes in prolific Garole sheep.Trop. Anim. Health Prod.42985–993. 10.1007/s11250-009-9518-1

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 55

  PruJ. K.TillyJ. L. (2001). Programmed cell death in the ovary: insights and future prospects using genetic technologies.Mol. Endocrinol.15845–853. 10.1210/mend.15.6.0646

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 56

  PryceJ.DaetwylerH. (2012). Designing dairy cattle breeding schemes under genomic selection: a review of international research.Anim. Prod. Sci.52107–114.

  • Google Scholar

• 57

  PurcellS.NealeB.Todd-BrownK.ThomasL.FerreiraM. A.BenderD.et al (2007). PLINK: a tool set for whole-genome association and population-based linkage analyses.Am. J. Hum. Genet.81559–575. 10.1086/519795

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 58

  Ruiz-CortésZ. T. (2012). “Gonadal sex steroids: production, action and interactions in mammals,” in Steroids-From Physiology to Clinical Medicine, ed.OstojicS. M. (London: IntechOpen).

  • Google Scholar

• 59

  SaitouN.NeiM. (1987). The neighbor-joining method: a new method for reconstructing phylogenetic trees.Mol. Biol. Evol.4406–425.

  • Pubmed Abstract
  • Google Scholar

• 60

  ShiY.WangS.BaiS.HuangL.HouY. (2015). Postnatal ovarian development and its relationship with steroid hormone receptors in JiNing Grey goats.Anim. Reprod. Sci.15439–47. 10.1016/j.anireprosci.2015.01.001

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 61

  StrangerB. E.StahlE. A.RajT. (2011). Progress and promise of genome-wide association studies for human complex trait genetics.Genetics187367–383. 10.1534/genetics.110.120907

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 62

  SuF.GuoX.WangY.WangY.CaoG.JiangY. (2018). Genome-wide analysis on the landscape of transcriptomes and their relationship with DNA methylomes in the hypothalamus reveals genes related to sexual precocity in Jining gray goats.Front. Endocrinol.9:501. 10.3389/fendo.2018.00501

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 63

  TajimaF. (1989). Statistical method for testing the neutral mutation hypothesis by DNA polymorphism.Genetics123585–595.

  • Pubmed Abstract
  • Google Scholar

• 64

  TalebiR.AhmadiA.AfrazF.SarryJ.WoloszynF.FabreS. (2018). Detection of single nucleotide polymorphisms at major prolificacy genes in the Mehraban sheep and association with litter size.Ann. Anim. Sci.18685–698. 10.2478/aoas-2018-0014

  • CrossRef
  • Google Scholar

• 65

  ThomasF. H.IsmailR. S.JiangJ.-Y.VanderhydenB. C. (2008). Kit ligand 2 promotes murine oocyte growth in vitro.Biol. Reprod.78167–175. 10.1095/biolreprod.106.058529

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 66

  TsukamotoS.MizutaT.FujimotoM.OhteS.OsawaK.MiyamotoA.et al (2014). Smad9 is a new type of transcriptional regulator in bone morphogenetic protein signaling.Sci. Rep.4:7596. 10.1038/srep07596

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 67

  TurnerS. D. (2014). qqman: an R package for visualizing GWAS results using Q-Q and manhattan plots.Biorxiv [Preprint]

  • Google Scholar

• 68

  VågeD. I.HusdalM.KentM. P.KlemetsdalG.BomanI. A. (2013). A missense mutation in growth differentiation factor 9 (GDF9) is strongly associated with litter size in sheep.BMC Genet.14:1. 10.1186/1471-2156-14-1

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 69

  VanRadenP.Van TassellC.WiggansG.SonstegardT.SchnabelR.TaylorJ.et al (2009). Invited review: reliability of genomic predictions for North American Holstein bulls.J. Dairy Sci.9216–24. 10.3168/jds.2008-1514

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 70

  VittiJ. J.GrossmanS. R.SabetiP. C. (2013). Detecting natural selection in genomic data.Annu. Rev. Genet.4797–120. 10.1146/annurev-genet-111212-133526

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 71

  WangK.LiM.HakonarsonH. (2010). ANNOVAR: functional annotation of genetic variants from high-throughput sequencing data.Nucleic Acids Res.38:e164. 10.1093/nar/gkq603

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 72

  XuL.ZhaoG.YangL.ZhuB.ChenY.ZhangL.et al (2019). Genomic patterns of homozygosity in Chinese local cattle.Sci. Rep.9:16977. 10.1038/s41598-019-53274-3

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 73

  XuS.-S.GaoL.XieX.-L.RenY.-L.ShenZ.-Q.WangF.et al (2018). Genome-wide association analyses highlight the potential for different genetic mechanisms for litter size among sheep breeds.Front. Genet.9:118. 10.3389/fgene.2018.00118

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 74

  YangJ.LeeS. H.GoddardM. E.VisscherP. M. (2011). GCTA: a tool for genome-wide complex trait analysis.Am. J. Hum. Genet.8876–82. 10.1016/j.ajhg.2010.11.011

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 75

  YangW.YanH.WangK.CuiY.ZhouT.XuH.et al (2019). Goat PDGFRB: unique mRNA expression profile in gonad and significant association between genetic variation and litter size.R. Soc. Open Sci.6:180805. 10.1098/rsos.180805

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 76

  YapD. B.WalkerD. C.PrenticeL. M.MckinneyS.TurashviliG.Mooslehner-AllenK.et al (2011). Mll5 is required for normal spermatogenesis.PLoS One6:e27127. 10.1371/journal.pone.0027127

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 77

  YuD.ChenF.ZhangL.WangH.ChenJ.ZhangZ.et al (2017). Smad9 is a key player of follicular selection in goose via keeping the balance of LHR transcription.bioRxiv [Preprint]

  • Google Scholar

• 78

  YuY.OuyangY.YaoW. (2017). shinyCircos: an R/Shiny application for interactive creation of Circos plot.Bioinformatics341229–1231. 10.1093/bioinformatics/btx763

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 79

  ZhangC.DongS.-S.XuJ.-Y.HeW.-M.YangT.-L. (2018). PopLDdecay: a fast and effective tool for linkage disequilibrium decay analysis based on variant call format files.Bioinformatics351786–1788. 10.1093/bioinformatics/bty875

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 80

  ZhangR.-Q.LaiF.-N.WangJ.-J.ZhaiH.-L.ZhaoY.SunY.-J.et al (2018). Analysis of the SNP loci around transcription start sites related to goat fecundity trait base on whole genome resequencing.Gene6431–6. 10.1016/j.gene.2017.12.002

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 81

  ZhangR.-Q.WangJ.-J.ZhangT.ZhaiH.-L.ShenW. (2019). Copy-number variation in goat genome sequence: a comparative analysis of the different litter size trait groups.Gene69640–46. 10.1016/j.gene.2019.02.027

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 82

  ZhangX.YanH.WangK.ZhouT.ChenM.ZhuH.et al (2018). Goat CTNNB1: mRNA expression profile of alternative splicing in testis and association analysis with litter size.Gene679297–304. 10.1016/j.gene.2018.08.061

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 83

  ZhaoZ.ParkC.McdevittM. A.Glidewell-KenneyC.ChambonP.WeissJ.et al (2009). p21-Activated kinase mediates rapid estradiol-negative feedback actions in the reproductive axis.Proc. Natl. Acad. Sci. U.S.A.1067221–7226. 10.1073/pnas.0812597106

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 84

  ZhouY.ZhouB.PacheL.ChangM.KhodabakhshiA. H.TanaseichukO.et al (2019). Metascape provides a biologist-oriented resource for the analysis of systems-level datasets.Nat. Commun.10:1523. 10.1038/s41467-019-09234-6

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 85

  ZhuH.ZhangY.BaiY.YangH.YanH.LiuJ.et al (2019). Relationship between SNPs of POU1F1 gene and litter size and growth traits in shaanbei white cashmere goats.Animals9:114. 10.3390/ani9030114

  • Pubmed Abstract
  • CrossRef
  • Google Scholar