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

Di Zhou1, Qingmeng Long1*, Rong Yang1, Xiaoshan Tan1, Jun Li1, Mingyan Tang1, Zhonghai Zhao3, Ye Ao2, Zhinan Zhou1 and Changxue Chen1

1Guizhou Testing Center for Livestock and Poultry Germplasm, Guiyang 550018, China

2Zunyi Animal Husbandry and Fishery Station, Zunyi, Guizhou, China, 563121

ABSTRACT

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.

Article Information

Received 18 December 2021

Revised 12 February 2022

Accepted 03 March 2022

Available online 22 September 2022

(early access)

Published 13 October 2023

Authors’ Contribution

DZ designed the project, and drafted the manuscript. QL provided experimental ideas. JL, XT and MT performed the experiment. RY, ZZ, YA and CC analyzed the transcriptome data.

Key words

Transcriptome, Black goat, Lambing

DOI: https://dx.doi.org/10.17582/journal.pjz/20211218031227

* Corresponding author: L_qmeng@Sina.com

030-9923/2023/0006-2761 $ 9.00/0

Copyright 2023 by the authors. Licensee Zoological Society of Pakistan.

This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

Introduction

The 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 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 (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 up-regulated 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.

The qRT-PCR amplifications were conducted in a total volume of 20 μL that included 10 μL 2 × UltraSYBR Mix, 0.5 μL each primer (10 nmol) and 8.5 μL ddH2O. The amplification was performed as follows: 95 ℃ for 3 min, 95 ℃ for 30 s, Tm for the 30s, and 72 ℃ 30s for 40 cycles.

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).

Results

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).

 

Table I. RT-PCR primers used in this study.

Gene name	Sequence (5’-3’)	Tm (℃)	Product (bp)
AMH	CTCTGCGTGAGCTGAGCGT TTGCCGTAGCGTGGGTTG	56.3	127
CYP1- 9A1	TGGACAGGTTGGAGGAGGT CCAGGAAGAGGATGTTAGAGG	58.3	99
FOXL2	CGCAGAAGCCCCCATACTC CGGCACCTTGATGAAGCACT	60	194
FSHR	ATCATGTTGGTGGGCTGGGTCT TGCGCTTGGCTATCTTGGTGTC	58.3	256
GATA4	CCAGAGCGAAGAAGGCAAACAAG GGGACTGACCACAGGGGAAGAGAA	59.4	140
INHBA	GGGTACGTGGAGATAGAGGAC TTGAAAGAGACGGATGGTGA	58.3	237
LHB	TTCACCACCAGCATCTGCG GAGACCATTGGGTCCACG	62.6	159
LHX9	AGACTCCGTCTACCACCTGA TGCAAGAGGGTCTCGAAGTG	55.5	126
ZP3	TGTGGTTAGGTTTGTGGTCGG TGGGTGTCATCTTCTCGGCG	49	302
CNOT1	TAAGGCTTGCAGAGGTAGGG TAAGTTGGCGGATTGAGGGA	56.3	240
DMRTA1	TCAGGAGGGGAAGAGAGTC TGTGAAGGGGAGAGAAAGC	55.3	385
GHRHR	TGACCGCCTCTCATTTCGCTACCAT CCAGCCACCAGAAGACCCTCCTTGT	58.4	120
Hox-A11	AGCGTGGTCCCTGCTCCTCTAAC GGGCTCAATGGCGTATTCTCTGA	60.3	194
Oviductin	GGGAAGTACACCAAGCAAGC TGTAACCGAAGCTGATGGCA	62.6	157
PCYT1B	ACTTGCTGGTAGGAGTTTG TTCTGTTCTCTGTGTTGGA	53.6	275
SLC34A2	TGATGAGTCGGTCCAAAACA GGATGACGCCAACAATAGC	61.4	236
WDR19	ATTGCATACCTCACCTCCCT GATTATTCATTCCCACAGCC	60.3	139
Wnt-7A	GCCTGGACGAGTGTCAGTTTCA ATAGCGGATGTCGGCAGAGCA	56.3	269
GAPDH	GCAAGTTCCACGGCACAG TCAGCACCAGCATCACCC	58.3	118

 

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).

 

Table II. Summary of read numbers from tissue samples used for this study.

Sample	Total clean bases (Gb)	Clean reads Q20(%)	Clean reads Q30(%)	Clean reads Ratio(%)	Total clean reads	Total mapping ratio (%)
ML_ovary_1	8.41	96.06	87.12	93.65	56049558	90.87
ML_ovary_2	9.04	96.13	87.31	92.58	60272540	90.95
ML_ovary_3	7.47	96.29	87.8	92.19	49820198	91.22
ML_oviduct_1	8.42	96.03	86.98	94.12	58154278	91.73
ML_oviduct_2	8.24	96.05	87.02	93.05	54907864	91.73
ML_oviduct_3	8.57	96.2	87.44	92.44	57116654	91.23
ML_uterus_1	7.37	95.9	86.77	92.84	49107842	90.56
ML_uterus_2	8.9	95.87	86.77	92.22	59325270	90.61
ML_uterus_3	8.5	96.34	87.83	92.71	53662236	90.90
SL_ovary_1	8.87	96.11	87.2	93.53	59157636	91.38
SL_ovary_2	8.06	96.04	87.09	93.49	53746050	91.31
SL_ovary_3	8.84	95.92	86.94	93.24	58904694	90.06
SL_oviduct_1	8.14	96.02	87	91.13	54263774	90.82
SL_oviduct_2	8.27	96.19	87.54	91.75	55135456	90.59
SL_oviduct_3	7.39	96.19	87.56	93.79	49234294	91.22
SL_uterus_1	8.36	96.78	86.52	93.78	55750482	91.04
SL_uterus_2	7.14	96.02	87.02	93.71	47612874	91.05
SL_uterus_3	11.16	95.94	86.83	94.05	74382518	91.18

 

SL, single lamb; ML, muliple lamb.

 

 

Figure 1 shows analysis of gene expression. An analysis of expression between groups indicated the presence of 29522 expressed genes in ovarian tissue (Fig. 2A), 28741 in fallopian tubes (Fig. 2B) and 28948 in uterine tissues (Fig. 2C).

 

 

Table III. Analysis of transcriptome results.

Sample	Total gene	Known gene	Novel gene	Total transcript	Known transcript	Novel transcript
ML_ovary_1	26220	25048	1172	39163	20399	18764
ML_ovary_2	27557	26339	1218	41773	22286	19487
ML_ovary_3	26963	25832	1131	40360	21701	18659
ML_oviduct_1	26455	25243	1212	40425	20926	19499
ML_oviduct_2	26050	24840	1210	39867	20401	19466
ML_oviduct_3	27040	25806	1234	41409	21759	19650
ML_uterus_1	26474	25304	1170	39918	20963	18955
ML_uterus_2	27409	26240	1169	41227	22041	19186
ML_uterus_3	26335	25136	1199	39594	20595	18999
SL_ovary_1	27749	26495	1254	42404	22461	19943
SL_ovary_2	28152	26906	1246	42936	23042	19894
SL_ovary_3	25264	24158	1106	37171	19372	17799
SL_oviduct_1	27113	25897	1216	41545	22072	19473
SL_oviduct_2	26657	25462	1195	40675	21271	19404
SL_oviduct_3	27203	26013	1190	41345	22001	19344
SL_uterus_1	27140	25937	1203	41220	21758	19462
SL_uterus_2	26923	25735	1188	40782	21472	19310
SL_uterus_3	27638	26420	1218	41976	22303	19673

 

SL, single lamb; ML, muliple lamb.

 

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 up-regulated, 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 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 up-regulated 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 (Zhang et al., 2021), sheep (Miao and Qin, 2015) and goats (Hou et al., 2015). For example, the gene RNAoar_circ_0012110 participates in the lambing traits of sheep by altering the activity of gonadotropin-releasing hormone (GnRH) as well as affecting the expression of key genes for this pathway (Zhang et al., 2019). 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 (Li et al., 2021). 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). Forked head L2 (FOXL2) is a member of the FOXL class 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, XM018049897.1, NM001285636.1, XM018051833.1, XM013963334.2, XM018061996.1, XM018061996.1, XM018062755.1, XM018060812.1, XM018041042.1, XM005692051.3, XM13965867.2, XM005679238.3, XM018047092.1, XM018044011.1, XM005681484.3, XM018055086.1, XM018049377.1, XM018066672.1 or available from the authors.

Funding

The project was supported by the Central Guidance to Local Funds, Grant (QKZY-2018-4015); the Core Program of Science and Technology Department of Guizhou Province (QKH[2019]2279)

Disclosure statement

The authors have declared no conflict of interest.

References

Ahlawat, S., Sharmaa, R., Roy, M., Mandakmale, S., Prakash, V. and Tantia, M.S., 2016. Genotyping of novel SNPs in BMPR1B, BMP15, and GDF9 genes for association with prolificacy in seven Indian goat breeds. Anim. Biotechnol., 27: 199-207. https://doi.org/10.1080/10495398.2016.1167706

Almeida, F.R.C.L., Costermans, N.G.J., Soede, N.M., Bunschoten, A., Keijer, J., Kemp, B. and Teerds, K.J., 2018. Presence of anti-Müllerian hormone (AMH) during follicular development in the porcine ovary. PLoS One, 13: 1–18. https://doi.org/10.1371/journal.pone.0197894

Amirpour, N.H., Khansefid, M., Mahmoud, G.G., Zhou, H. and Hickford, J.G.H., 2020. Identification of polymorphisms in the oocyte-derived growth differentiation growth factor 9 (GDF9) gene associated with litter size in New Zealand sheep (Ovis aries) breeds. Reprod. Domest. Anim., 55: 1585–1591. https://doi.org/10.1111/rda.13813

An, X.P., Hou, J.X., Lei, Y.N., Gao, T.Y., Song, Y.X., Wang, J.G. and Cao, B.Y., 2015. Two mutations in the 5′-flanking region of the KITLG gene are associated with litter size of dairy goats. Anim. Genet., 46: 308–311. https://doi.org/10.1111/age.12277

Bai, H., Sakurai, T., Godkin, J.D. and Imakawa, K., 2014. Expression and in situ localization of GATA4, 5 and 6mRNAs in ovine conceptuses and uterine endometria during the peri-implantation period. Anim. Sci. J., 85: 388–394. https://doi.org/10.1111/asj.12156

Blanca, A., Verónica, M., Manuel, A., Alfonso, G.A., and Dimitrios, R.M.J.M., 2018. Effects of recombinant OVGP1 protein on in vitro bovine embryo development. J. Reprod. Dev., 64: 433-443. https://doi.org/10.1262/jrd.2018-058

Cadoret, V., Frapsauce, C., Jarrier, P., Maillard, V., Bonnet, A., Locatelli, Y., Royère, D., Monniaux, D., Guérif, F. and Monget, P., 2017. Molecular evidence that follicle development is accelerated in vitro compared to in vivo. Reproduction, 153: 493–508. https://doi.org/10.1530/REP-16-0627

Chantepie, L., Bodin, L., Sarry, J., Woloszyn, F., Plisson-Petit, F., Ruesche, J., Drouilhet, L. and Fabre, S. 2020. Genome-wide identification of a regulatory mutation in BMP15 controlling prolificacy in sheep. Front. Genet., 11: 585. https://doi.org/10.3389/fgene.2020.00585

Chong, Y., Huang, H., Liu, G., Jiang, X. and Rong, W., 2018. A single nucleotide polymorphism in the zona pellucida 3 gene is associated with the first parity litter size in Hu sheep. Anim. Reprod. Sci., 193: 26–32. https://doi.org/10.1016/j.anireprosci.2018.03.028

Cleys, E.R., Halleran, J.L., Enriquez, V.A., Da Silveira, J.C., West, R.C., Winger, Q.A., Anthony, R.V., Bruemmer, J.E., Clay, C.M. and Bouma, G.J., 2015. Androgen receptor and histone lysine demethylases in ovine placenta e0117472. PLoS One, 10: 1–19. https://doi.org/10.1371/journal.pone.0117472

De Castro, F., Ruiz, R., Montoro, L., Pérez-Hernández, D., Sánchez-Casas Padilla, E., Real, L.M. and Ruiz, A., 2003. Role of follicle-stimulating hormone receptor Ser680Asn polymorphism in the efficacy of follicle-stimulating hormone. Fertil. Steril., 80: 571–576. https://doi.org/10.1016/S0015-0282(03)00795-7

Dutta, R., Laskar, S., Borah, P., Kalita, D., Das, B., Zaman, G., Barman, N.N. and Saikia, D.P., 2014. Polymorphism and nucleotide sequencing of BMPR1B gene in prolific Assam hill goat. Mol. Biol. Rep., 41: 3677–3681. https://doi.org/10.1007/s11033-014-3232-4

Eozenou, C., Lesage-padilla, A., Mauffré, V., Healey, G.D., Camous, S., Bolifraud, P., Giraud-delville, C., Vaiman, D., Shimizu, T., Miyamoto, A., Sheldon, I.M., Constant, F., Pannetier, M. and Sandra, O., 2020. FOXL2 is a progesterone target gene in the endometrium of ruminants. Int. J. mol. Sci., 21: 1478. https://doi.org/10.3390/ijms21041478

Green, M.R., and Sambrook, J., 2019. Isolation of poly (A)+ messenger RNA using magnetic oligo (dT) beads. Cold Spring Harbor Protocols, 2019, 2019(10): pdb. prot101733. https://doi.org/10.1101/pdb.prot101733

Gunawan, A., Sahadevan, S., Neuhoff, C., Große-Brinkhaus, C., Gad, A., Frieden, L., Tesfaye, D., Tholen, E., Looft, C., Uddin, M.J., Schellander, K. and Cinar, M.U., 2013. RNA deep sequencing reveals novel candidate genes and polymorphisms in boar testis and liver tissues with divergent androstenone levels. PLoS One, 8: e63259. https://doi.org/10.1371/journal.pone.0063259

Guo, J., Schally, A.V., Zarandi, M., Varga, J. and Leung, P.C.K., 2010. Antiproliferative effect of growth hormone-releasing hormone (GHRH) antagonist on ovarian cancer cells through the EGFR-Akt pathway. Reprod. Biol. Endocrinol., 8: 1–11. https://doi.org/10.1186/1477-7827-8-54

Hernández-Montiel, W., Collí-Dula, R.C., Ramón-Ugalde, J.P., Martínez-Núñez, M.A. and Zamora-Bustillos, R., 2019. RNA-seq transcriptome analysis in ovarian tissue of pelibuey breed to explore the regulation of prolificacy. Genes (Basel), 10: 1–13. https://doi.org/10.3390/genes10050358

Hou, J.X., An, X.P., Han, P., Peng, J.Y. and Cao, B.Y., 2015. Two missense mutations in exon 9 of caprine PRLR gene were Asso. with litter size. Anim. Genet., 46: 87–90. https://doi.org/10.1111/age.12223

Ito, K., Takahashi, A., Morita, M., Suzuki, T. and Yamamoto, T., 2011. The role of the CNOT1 subunit of the CCR4-NOT complex in mRNA deadenylation and cell viability. Protein Cell, 2: 755–763. https://doi.org/10.1007/s13238-011-1092-4

Knarston, I.M., Pachernegg, S., Robevska, G., Ghobrial, I., Er, P.X., Georges, E., Takasato, M., Combes, A.N., Jørgensen, A., Little, M.H., Sinclair, A.H. and Ayers, K.L., 2020. An in vitro differentiation protocol for human embryonic bipotential gonad and testis cell development. Stem. Cell Rep., 15: 1377–1391. https://doi.org/10.1016/j.stemcr.2020.10.009

Kwon, S.G., Hwang, J.H., Park, D.H., Kim, T.W., Kang, D.G., Kang, K.H., Kim, I.S., Park, H.C., Na, C.S., Ha, J. and Kim, C.W., 2016. Identification of differentially expressed genes associated with litter size in Berkshire pig placenta. PLoS One, 11: 1–12. https://doi.org/10.1371/journal.pone.0153311

La, Y., Tang, J., He, X., Di, R., Wang, X., Liu, Q., Zhang, L., Zhang, X., Zhang, J., Hu, W. and Chu, M., 2019. Identification and characterization of mRNAs and lncRNAs in the uterus of polytocous and monotocous Small Tail Han sheep (Ovis aries). Peer J., 2019: 1–21. https://doi.org/10.7717/peerj.6938

Lai, L., Shen, X., Liang, H., Deng, Y., Gong, Z. and Wei, S., 2018. Determine the role of FSH receptor binding inhibitor in regulating ovarian follicles development and expression of FSHR and ER α in mice. Biomed. Res. Int., 2018: Article ID 5032875. https://doi.org/10.1155/2018/5032875

Li, C., He, X., Zhang, Z., Ren, C. and Chu, M., 2021. Pineal gland transcriptomic profiling reveals the differential regulation of lncRNA and mRNA related to prolificacy in STH sheep with two FecB genotypes. BMC Genom. Data, 22: 1–17. https://doi.org/10.1186/s12863-020-00957-w

Mahdavi, M., Nanekarani, S. and Hosseini, S.D., 2014. Mutation in BMPR-IB gene is associated with litter size in Iranian Kalehkoohi sheep. Anim. Reprod. Sci., 147: 93–98. https://doi.org/10.1016/j.anireprosci.2014.04.003

Martinez-Royo, A., Jurado, J.J., Smulders, J.P., Martí, J.I., Alabart, J.L., Roche, A., Fantova, E., Bodin, L., Mulsant, P., Serrano, M., Folch, J. and Calvo, J.H., 2008. A deletion in the bone morphogenetic protein 15 gene causes sterility and increased prolificacy in Rasa Aragonesa sheep. Anim. Genet., 39: 294–297. https://doi.org/10.1111/j.1365-2052.2008.01707.x

Mecha, M., Yanguas-Casás, N., Feliú, A., Mestre, L., Carrillo-Salinas, F.J., Riecken, K., Gomez-Nicola, D. and Guaza, C., 2020. Involvement of Wnt7a in the role of M2c microglia in neural stem cell oligodendrogenesis. J. Neuroinflam., 17: 1–17. https://doi.org/10.1186/s12974-020-01734-3

Miao, X. and Qin, Q.L.X., 2015. Genome-wide transcriptome analysis of mRNAs and microRNAs in Dorset and Small Tail Han sheep to explore the regulation of fecundity. Mol. Cell. Endocrinol., 402: 32–42. https://doi.org/10.1016/j.mce.2014.12.023

Miao, X., Luo, Q., Zhao, H. and Qin, X., 2016. Co-expression analysis and identification of fecundity-related long non-coding RNAs in sheep ovaries. Sci. Rep., 6: 1–10. https://doi.org/10.1038/srep39398

Nicol, B., Rodriguez, K. and Yao, H.H.C., 2020. Aberrant and constitutive expression of FOXL2 impairs ovarian development and functions in mice. Biol. Reprod., 103: 966–977. https://doi.org/10.1093/biolre/ioaa146

Pokharel, K., Peippo, J., Honkatukia, M., Seppälä, A., Rautiainen, J., Ghanem, N., Hamama, T.M., Crowe, M.A., Andersson, M., Li, M.H. and Kantanen, J., 2018. Integrated ovarian mRNA and miRNA transcriptome profiling characterizes the genetic basis of prolificacy traits in sheep (Ovis aries). BMC Genom., 19: 1–17. https://doi.org/10.1186/s12864-017-4400-4

Rajani, M., George, Katherine, L., Rawls, A., Rawls, A., Viger, R. and Wilson-Rawls. 2015. Notch signaling represses GATA4-induced expression of genes involved in steroid biosynthesis. Reproduction, 150: 383–394. https://doi.org/10.1530/REP-15-0226

Rather, M.A., Bhat, I.A. and Sharma, R., 2016. Identification, cDNA cloning, and characterization of luteinizing hormone beta subunit (Lhb) gene in Catla catla. Anim. Biotechnol., 27: 148–156. https://doi.org/10.1080/10495398.2016.1140055

Ren, J., Sun, C., Chen, L., Hu, J., Huang, X., Liu, X. and Lu, L., 2018. Exploring differentially expressed key genes related to development of follicle by RNA-seq in Peking ducks (Anas platyrhynchos). PLoS One, 14: e0209061. https://doi.org/10.1101/483818

Sciences, V., 2017. Clinical study with investigation of polymorphism FSHR gene in twining goats in Iraq. Adv. Anim. Vet. Sci., 5: 192–195.

Thomas, N., Venkatachalapathy, T., Aravindakshan, T. and Raghavan, K.C., 2016. Molecular cloning, SNP detection and association analysis of 5’ flanking region of the goat IGF1 gene with prolificacy. Anim. Reprod. Sci., 167: 8–15. https://doi.org/10.1016/j.anireprosci.2016.01.016

Torres-Rovira, L., Succu, S., Pasciu, V., Manca, M.E., Gonzalez-Bulnes, A., Leoni, G.G., Pennino, M.G., Spezzigu, A., Gallus, M., Dattena, M., Monniaux, D., Naitana, S. and Berlinguer, F., 2016. Postnatal pituitary and follicular activation: A revisited hypothesis in a sheep model. Reproduction, 151: 215–225. https://doi.org/10.1530/REP-15-0316

Våge, D.I., Husdal, M., Kent, M.P., Klemetsdal, G. and Boman, I.A., 2013. A missense mutation in growth differentiation factor 9 (GDF9) is strongly associated with litter size in sheep. BMC Genet., 14: 1–8. https://doi.org/10.1186/1471-2156-14-1

Xu, S., Xia, W., Zohar, Y. and Gui, J.F., 2013. Zebrafish dmrta2 regulates the expression of cdkn2c in spermatogenesis in the adult testis. Biol. Reprod., 88: 1–12. https://doi.org/10.1095/biolreprod.112.105130

Yue, B., Yang, H., Wu, J., Wang, J., Ru, W., Cheng, J., Huang, Y., Lei, C., Lan, X. and Chen, H., 2020. Characterization and transcriptome analysis of exosomal and nonexosomal rnas in bovine adipocytes. Int. J. mol. Sci., 21: 1–17. https://doi.org/10.3390/ijms21239313

Zhang, G.M., Guo, Y.X., Cheng, C.Y., El-Samahy, M.A., Tong, R., Gao, X.X., Deng, K.P., Wang, F. and Lei, Z., 2020. Arginine infusion rescues ovarian follicular development in feed-restricted Hu sheep during the luteal phase. Theriogenology, 158: 75–83. https://doi.org/10.1016/j.theriogenology.2020.09.002

Zhang, X., Wang, R., Wang, T., Zhang, X., Dongye, M., Wang, D., Wang, J., Li, W., Wu, X., Lin, D. and Lin, H., 2021. The metabolic reprogramming of Frem2 mutant mice embryos in cryptophthalmos development. Front. Cell Dev. Biol., 8: 1–9. https://doi.org/10.3389/fcell.2020.625492

Zhang, Z., Tang, J., He, X., Zhu, M., Gan, S., Guo, X., Zhang, X., Zhang, J., Hu, W. and Chu, M., 2019. Comparative transcriptomics identify key hypothalamic circular rnas that participate in sheep (Ovis aries) reproduction. Animals, 9: 557. https://doi.org/10.3390/ani9080557

Zi, X.D., Hu, L., Lu, J.Y., Liu, S. and Zheng, Y.C., 2020. Comparison of the sequences and expression levels of genes related to follicular development and atresia between prolific and nonprolific goat breeds. Vet. Med. Sci., 6: 187–195. https://doi.org/10.1002/vms3.225