Screening of Candidate Genes and Tissue Expression Analysis of Lambing Traits for the Guizhou Black Goat

We used RNAseq analysis to identify key genes affecting reproductive traits for the Guizhou black goat. We utilized mRNA collected from uterine, ovarian and fallopian tube tissues of goats with the same pedigree and analyzed the results using bioinformatics and qRT-PCR. We found 8120 differentially expressed genes (DEGs) in ovarian tissues (3205 up- and 4915 down-regulated), 5255 DEGs in the fallopian tubes (1317 up- and 3888 down-regulated) and 5180 DEGs in uterine tissue (2597 up- and 2583 down-regulated). The biological processes associated with these DEGs were primarily cell differentiation, reproduction and biological regulation. The associated gene functions were related to organelle, cell membrane and cell structure in cell component processes and for molecular function, the genes were related to protein binding, catalytic activity and cell proliferation. Moreover, these DEGs were enriched in the categories of signal transduction, reproductive system, endocrine system and cancer-related pathways and linked the cGMP-PKG, cAMP, TGFβ and mTOR signaling pathways to reproductive performance. We verified these results using qRT-PCR for 9 up-regulated (AMH, CYP19A1, FOXL2, FSHR, GATA4, INHBA, LHB, LHX9 and ZP3) and 9 down-regulated (CNOT1. DMRTA1, GHRHR, Hox-All, Oviductin, PCYT1B, SLC34A2, WDR19 andWnt-7A) genes. The results were consistent with the RNAseq results and we more specifically linked AMH, GATA4, Oviductin and DMRTA1 to lambing traits of the Guizhou black goat.


INTRODUCTION
T he lambing trait is an important economic trait for goats and its genetics are an important basis for determining lamb numbers and improving goat production efficiency. Reproductive traits are currently assessed for candidate genes using the polymerase chain reaction (PCR) single strand confirmation polymorphism (SSCP) and restriction fragment linked polymorphism (RFLP) techniques. The relationships between different genotypes and reproductive traits are then analyzed to determine whether the candidate genes are directly involved in the trait and are therefore primary genes involved in the phenotype (Ahlawat et al., 2016;Amirpour Najafabadi et al., 2020). These methods have been used to screen candidate genes affecting goat lambing performance but the primary genes affecting production performance have not yet been confirmed (Dutta et al., 2014;An et al., 2015;Hou et al., 2015;Thomas et al., 2016).
The identification of lambing traits can also be analyzed using high-throughput sequencing technologies where transcriptiome, proteome or metabolome expression patterns of differentially expressed genes (DEGs) are used to link reproductive functions with molecular mechanism. For example, transcriptomics has been used for identification of related genes in pigs, sheep and goats and the primary genes affecting litter size in sheep have been successfully identified using this method. In particular, specific mRNA and miRNAs expressed in the sheep ovary have been linked to reproductive traits and 79 novel and relevant miRNAs have been identified (Pokharel et al., 2018). In addition, the FecB and BMPR-IB genes were linked to litter size in Kalehkoohi sheep and the FecX allele for the BMP15 mutation controlled the ovulation rate and reproductive ability. FecX homozygous ewes possessed strong reproductive ability and this was linked to decreased BMP15 promoter activity (Mahdavi et al., 2014;Chantepie et al., 2020). The sheep GDF9 gene is another candidate gene that linked ovulation rate and litter size (Våge et al., 2013). Polymorphisms of BMPR-IB have also been are associated with lambing traits in goats (Hou O n l i n e F i r s t A r t i c l e et al., 2015) as have single nucleotide polymorphisms (SNPs) in the 5 'regulatory region of the IGF1 gene (Thomas et al., 2016). Additionally, the RNAseq technique has been applied to the placenta of Berkshire pigs where the IL-6 and LIPG genes were identified as important in placental nutrition supply during pregnancy (Kwon et al., 2016). Transcriptome studies have also identified lncRNA targets of the TGFβ and oxytocin genes in sheep with FecB (HanBB and Han++) genotypes that enhanced sheep reproduction (Miao et al., 2016). The Guizhou black goat is primarily raised in Bijie, Guizhou, China and is a large animal with a pure black coat and a high meat quality, but suffers from low fecundity with an average lambing level of 152%. In efforts to improve the fecundity of Guizhou Black goats, we screened goats of similar pedigrees that could link litter traits using RNAseq analysis to identify DEGs and performed a functional and pathway enrichment analysis. The qRT-PCR method was used to confirm the differentially expressed genes in tissues and our result provide an experimental foundation for further research into Guizhou black goat lambing traits.

MATERIALS AND METHODS
All experiments followed the relevant guidelines and regulations formulated by the Ministry of Agriculture of the people's Republic of China. The animal handling procedures conformed with the China Animal Welfare guidelines and have been approved by the Guizhou Provincial Animal Protection and Use Commission (Guiyang, China). Six Guizhou black goats with lambing records and health conditions were sacrificed for this study and their uterus, ovary and fallopian tubes were collected for transcriptome analysis and follow-up functional verification. In this study, 6 Guizhou black goats with lambing records and good health status were selected as experimental animals in the same feeding environment. The twin lambs were 3.5 years old, and the birth was recorded as the second child. After sacrifice, the uterus, ovaries and fallopian tubes were collected for transcriptome analysis and subsequent functional verification.

Reagents
Trizol reagent was obtained from Invitrogen (Carlsbad, CA, USA). Chloroform and absolute ethanol were of reagent grade and obtained from local suppliers. The following were obtained from commercial suppliers; reverse transcription kit (company and location) 2 × Taq PCR Master Mix (company and location), UltraSYBR Mix (CWBIO, give location), RNeasy RNA purification kit (Qiagen, Valencia, CA, USA), Agilent RNA 6000 Nano Kit (Santa Clara, CA, USA), Oligo dT magnetic bead mRNA purification kit (company and location), polynucleotide kinase (company and location) and DL2 molecular size markers (Takara, Dalian, China).

Sequencing library construction
RNA was extracted from was extracted from the uterus, fallopian tubes and ovary tissues using Trizol and further purified using the RNAeasy kit using the recommendations of the manufacturers. RNA quantity and quality were assessed using an Agilent 2100 Bioanalyzer with an RNA 6000 Nano Kit as suggested by the manufacturer. Poly A mRNA was isolated from total RNA using oligo (dT) coupled to magnetic beads that were then fragmented using ion interruption (Green et al., 2019). The short fragments were then used as templates for first strand cDNA synthesis using 6 bp random primers for reverse transcription and second strand cDNA was synthesized using Taq DNA polymerase using standard conditions. The ds DNA was phosphorylated at the 5 end using polynucleotide kinase and AT tailed at the 3 end, and then connecting a bubble-shaped linker with a protruding "T" at the 3' end. The linked product was used to generate ssDNA using xxxx polymerase and bridge primers using rolling circle synthesis as previously reported. The resulting ssDNA circular library was used for high throughput sequence analysis (see below).

Illumina sequence processing
The IlluminaHiSeq sequencing platform was used for DNA sequence analysis at a commercial company (Shenzhen Huada Genome, give location). The original offline sequence was processed using SOAPnuke software (https://github.com/BGI-flexlab/SOAPnuke) and the rawreads were processed by removing joining sequence information and sequences < 5bp. For double-end sequencing, if one end was contaminated by the connector, the read was removed as were low-quality reads whose quality value was < 10 and accounted for >20% of the total base number. Reads containing > 5% N were also eliminated.

Genome alignment and expression analysis
The clean reads were then mapped to the goat Capra hircus reference genome (https://www.ncbi.nlm.nih.gov/ assembly/GCF_001704415.1) using HISAT http:// www. ccb.jhu.edu/software/hisat/index.shtml. In particular, the RSEM module was used to calculate the expression level of each gene to obtain the fragments per kilobase of transcript per million mapper reads (FPKM) value and represented the level of gene expression.

Identification of differentially expressed genes
The FPKM value was analyzed using DEGseq

O n l i n e F i r s t A r t i c l e
(https://bioconductor.org/packages/release/bioc/html/ DEGseq.html) and the P-value was calculated according to the normal distribution and RSEM was used to compare differences (Benjamini, Y. and Hochberg, Y; Storey, J. and Tibshirani, R). The P and Q values were correlated using the FDR value calculated using the Possion distribution function in the program and genes with FDR ≤ 0.001 and Q ≤ 0.001 were deemed as significantly and differentially expressed.

GO enrichment and KEGG analysis of DEGs
The Gene Ontology database (GO: blue.Coexpression. network; KEGG: Modules_gene_bluePathway) was used for gene ontology (GO) and KEGG Pathway enrichment using the phyper function to calculate DEG numbers per term. The pathway was significantly enriched if the following conditions were met: FDR ≤ 0.01 and Q ≤ 0.001.

DEG verification
The DEG analyses were verified by selecting 9 upregulated and 9 down-regulated genes for verification using quantitative real-time PCR (qRT-PCR). PCR primers were designed using Primer-BLAST (https:// www.ncbi.nlm.nih.gov/tools/primer-blast/) (Table I) and the amplification reactions were performed in triplicate and the goat GAPDH (XM_005680968.3) gene was used as an internal reference. The relative expression level of the tested sample was obtained by the 2 -ΔΔCt method.

Statistical analysis
The qRT-PCR data were calculated and analyzed using the 2 -ΔΔCt method and the one-way ANOVA of the test data was analyzed using SPSS20.0 software (IBM, Chicago, Ill, USA).

Transcriptome sequencing
Over 18 tissue samples from the uterus, fallopian tubes and ovary of 3 animals generated an average of 8.38Gb per sample that were 91.03% matched to the reference genome. The uniform comparison level between sample data were 100% (Table II).
These sequencing results indicated that the average total transcript for each tissue for all 3 animals was 40415 and the unknown transcript level was 19185 with 26722 known and 1190 unknown genes. For each animal, the average total number of known transcripts was 41117 and unknown was 19366 and these groups possessed 27093 and 1201 unknown genes on average, respectively (Table  III).

Gene expression analysis
The overall distribution of expression trends per sample was approximately the same but did not completely overlap. This indicated the presence of differences between samples (Fig. 1).

DEG analysis
In this study we used the DEGseq method was used to analyze gene expression differences and genes were included in the analysis if logFoldChange was ≥ 2 and Q < 0. 001. We found 8120 DEGs (3205 up-regulated, 4915 down-regulated) in ovarian tissue, 5255 DEGs (1317 upregulated, 3888 down-regulated) in the fallopian tube and 5180 DEGs (2597 up-regulated, 2583 down-regulated) in uterine tissue (Fig. 3).

GO functional enrichment analysis of DEGs
Gene ontology functional items that were significantly enriched in the DEGs were corrected using the FDR and were regarded as significantly enriched if | logFoldChange | ≥ 2 and Q ≤ 0.001. Within these criteria we identified 2858, 1855 and 1836 genes that were enriched in the ovary, fallopian tube and uterus, respectively. The biological process grouping primarily consisted of cell differentiation, reproductive process and biological regulation; the cellular components included the primary groups structural component of organelles, cell membranes and cells. The molecular function group consisted of items relating to protein binding, catalytic activity and cell proliferation (Fig. 4A-C).

KEGG enrichment analysis of DEGs
All transcripts from ovary, uterus and fallopian tubes of single and multi-lamb groups were annotated to the KEGG database, and the number of differentially expressed genes at different levels were counted when the condition | logFoldChange | ≥ 2 and Q ≤ 0.001 was met. The metabolic and signaling pathways were then determined. We found that the ovary DEGs could be annotated to 334

D. Zhou et al.
KEGG pathways and included 4469 DEGs (Fig. 5A). The fallopian tube genes were annotated to 331 KEGG pathways and included 2789 DEGs (Fig. 5B). The transcripts from the uterus samples were annotated to 332 KEGG pathways and contained 2893 DEGs (Fig. 5C). These DEGs were primarily grouped in the signal transduction, reproductive system, endocrine system and cancer-related pathways. These included cGMP-PKG, cAMP, TGF-β and mTOR signaling pathways and their regulatory mechanisms await further investigation. However, we found a 1:1 correspondence of the qRT-PCR and RNAseq results and the levels of regulation were consistent (Fig. 6A-F). These results were verified using qRT-PCR of 9 upregulated and 9 down-regulated genes.

DISCUSSION
The RNAseq method of gene expression analysis has been, used for screening related genes involved in reproductive traits in pigs (Gunawan et al., 2013), cattle (Yue et al., 2020), mice , sheep (Miao and Qin, 2015) and goats . For example, the gene RNAoar_circ_0012110 participates in the lambing traits of sheep by altering the activity of gonadotropinreleasing hormone (GnRH) as well as affecting the expression of key genes for this pathway . The differentially expressed lncRNA in the pineal gland affects sheep yields by regulating steroid production and hormone activity via regulation of the target gene XLOC_466330 . The adenylate-uridylate (AU)-rich elements (AREs) in the 3' untranslated regions of transcripts in the placenta responding to androgen signaling identified the recruitment of the AR-KDM complex that regulated placental VEGFA expression and trophoblast differentiation and placental development (Cleys et al., 2015). Additionally, transcriptome sequencing analysis determined that the CA5A and FEM1B genes were involved in the final stage of follicular development (Hernández-Montiel et al., 2019). Follicular development was also regulated by EPOX in ducks during oviposition (Ren et al., 2018). Transcriptome sequencing technology can serve as an effective technology to mine unknown but related genes. At present, sheep proliferation genes have been successfully screened and applied to sheep breeding protocols using this approach and have identified the involvement of GDF9, FecX, FecB and BMPRIB (Martinez-Royo et al., 2008). The primary genes or effective molecular markers that determine lambing traits have not been identified for goats. Therefore, in this study we used the Guizhou black goat with multiple lambs as the research object and single lambs as the control group and applied the RNAseq methodology using uterus, ovary and fallopian tube tissues. We found the presence of both up-regulated and down-regulated genes involved in cell proliferation, differentiation, protein binding and signal transduction in the uterus, ovary and fallopian tubes. The genes related to lambing were screened from the sequencing results and verified by transcriptome data and the qRT-PCR results were consistent with that of the transcriptome. We could identify specific reproductive functions for the DEGs we identified in this study. For instance, our group of up-regulated genes include the anti-Mullerian hormone gene (AMH) expressed by developing follicular granulosa cells and its expression is highest in mature ovaries. AMH plays a regulatory role in follicular development by stimulating follicular hormone to threshold levels (Ren et al., 2018) and can promote FSH secretion and improve the conception rate of sheep (Almeida et al., 2018). The aromatase (CYP19A1) gene is crucial for follicular development and increases steroid production and is abundantly expressed in mature follicles (Torres-Rovira et al., 2016;Cadoret et al., 2017;Zhang et al., 2020). of transcription factors and is a key element for ovarian differentiation and function (Eozenou et al., 2020). It is also closely related to follicular development and granulosa cell proliferation and differentiation and reduced expression levels significantly reduce steroid production (Nicol et al., 2020). Another gene (FSHR) plays an important role in the regulation of steroid production and follicular proliferation and its deficiency inhibits follicular development (Lai et al., 2018). FSHR has been linked to the twin lamb traits of goats (Sciences, 2017) and regulates steroid production and follicular proliferation during ovarian maturation, animals with insufficient FSHR are likely to be infertile (De Castro et al., 2003). The zinc finger transcription factor 4 (GATA4) exists in uterine epithelial cells and is specifically expressed in ovarian granulosa cells and can positively regulate steroid production (Bai et al., 2014;Rajani et al., 2015). Luteinizing hormone β subunit (LHB) is primarily expressed in the pituitary and ovary but is also involved in follicular development. Interestingly, LHB is a lncRNA target and regulates reproduction and lies in the signal transduction pathway and biological process units we identified (Rather et al., 2016;La et al., 2019;Zi et al., 2020). LHX9 is a transcription factor that regulates a variety of developmental processes and has a protective effect on embryos and also functions as a gonadal marker. The absence of LHX9 expression inhibits gonadal development (Knarston et al., 2020). The zona pellucida protein 3 (ZP3) gene is the primary sperm receptor and acrosome reaction inducer and plays an important role for reproductive traits. SNP (rs401271989) in its coding region affects the first litter size of the goat (Chong et al., 2018). These up-regulated genes are involved in embryonic development, hormone regulation and follicular development and are linked to key reproductive traits. Therefore, our transcriptome sequencing results are reliable and the enhanced expression of these up-regulated genes may have an important effect on the lambing traits of the Guizhou black goat.
We also identified genes that were down-regulated such as CNOT1 that has not been previously linked to reproduction but its gene product plays an essential role in mRNA turnover in the CCR4-NOT deadenylase complex (Ito et al., 2011). DMRTA1 plays an independent role in gonadal development and its over-expression during embryogenesis leads to serious underdevelopment and high expression in the ovary can result in polyfollicular syndrome (Xu et al., 2013). The GHRHR encodes the receptor gene of auxin-releasing hormone (GHRH) and the latter reduces cell viability by inhibiting GHRHR expression (Guo et al., 2010). Tubal protein (Oviductin) is engaged during fertilization and the early stages of embryonic development and its expression does not alter the rate of embryo development but can enhance embryo quality. (Blanca et al., 2018). Wnt-7A is a cytokine released by microglia that controls neural stem cell survival, proliferation, migration and differentiation (Mecha et al., 2020).
The down-regulated genes we identified inhibit embryonic development and cell viability implying that they could affect goat lambing traits. The functions of Hox-All, PCYT1B, SLC34A2 and WDR19 genes in reproductive performance are unknown and need to be further studied.
The integrity of our RNAseq results were also verified using qRT-PCR analysis. The expression of up-regulated genes in the multiple lamb groups were significantly higher than for the single-lamb groups while the down-regulated genes in the single-lamb group were significantly higher than that in multi-lamb group (P < 0.05). In general, our group of up-regulated genes were linked to enhanced goat litter characteristics whereas the down-regulated genes are inhibitory when over-expressed. These results provide basic experimental data for further research into the molecular regulation of reproductive traits.

CONCLUSIONS
In this study we used transcriptome sequencing and identified genes that were up-and down regulated in ovarian, fallopian tube and uterine tissues. We identified the genes that were differentially expressed and could be linked to reproductive processes either directly or indirectly. The gene functions could be classified as involved in signal transduction, reproductive and endocrine systems and cancer related pathways. The reliability of our results was verified using qRTPCR for 9 up-regulated genes (AMH, CYP19A1, FOXL2, FSHR, GATA4, INHBA, LHB, LHX9, ZP3) and 9 down-regulated genes (CNOT1, DMRTA1, GHRHR, Hox-All, Oviductin, PCYT1B, SLC34A2, WDR19 and Wnt-7A). The up-regulated genes we identified were linked to follicular growth and development, ovarian granulosa cell proliferation and gonadal development and these functions were closely related to multi-lamb reproductive traits. The down-regulated genes have been shown in previous studies to have adverse effects on embryogenesis and cell viability when over-expressed. Our results lay a foundation for the follow-up study of the molecular regulation of genes involved in reproductive traits.

Data availability statement
The data that support the findings of this study are openly available at https://www.ncbi.nlm.nih.gov/ nucleotide/ with the accession numbers XM018050766.1,