Identification of differentially expressed long noncoding RNAs in the ovarian tissue of Shal and Sangsari ewes using RNA-seq.

BACKGROUND: The ovary has an important role in reproductive function. Animal reproduction is dominated by numerous coding genes and noncoding elements. Although long noncoding RNAs (LncRNAs) are important in biological activity, little is known about their role in the ovary and fertility.
METHODS: Three adult Shal ewes and three adult Sangsari ewes were used in this investigation. LncRNAs in ovarian tissue from two breeds were identified using bioinformatics analyses, and then target genes of LncRNAs were discovered. Target genes were annotated using the DAVID database, and their interactions were examined using the STRING database and Cytoscape software. The expression levels of seven LncRNAs with their target genes were assessed by real-time PCR to confirm the RNA-seq.
RESULTS: Among all the identified LncRNAs, 124 LncRNAs were detected with different expression levels between the two breeds (FDR < 0.05). According to the DAVID database, target genes were discovered to be engaged in one biological process, one cellular component, and 21 KEGG pathways (FDR < 0.05). The PES1, RPS9, EF-1, Plectin, SURF6, CYC1, PRKACA MAPK1, ITGB2 and BRD2 genes were some of the most crucial target genes (hub genes) in the ovary.
CONCLUSION: These results could pave the way for future efforts to address sheep prolificacy barriers.

INTRODUCTION

Sheep are kept on farms to generate meat, milk, wool and fur, all of which are valuable human goods (Montossi et al., 2013). Fertility factors that influence sheep production include ovulation rate, litter size, total lambs born, age at first lambing, stillbirth, and age at first maturity (Chen et al., 2015). Reproductive traits differ not only between sheep breeds but also within the same flock. An important economic criterion is the classification of sheep with the highest prolificacy and fecundity. While phenotypic selection can be used to identify breeding stocks with a higher likelihood of increasing the reproductive rate, identifying relevant genes connected to prolificacy and fecundity can help with effective breeding (Gholizadeh et al., 2014). In terms of prolificacy, there are various differences among Iranian sheep breeds (Eghbalsaied et al., 2012). Most sheep breeds have one lamb, but only a few have twins or multiples (Miao et al., 2016a), and this trait is influenced by both environmental and genetic factors (Drouilhet et al., 2009). Multiple major genes, such as GDF9, BMPR1B, B4GALNT2 and FecX2, influence ovulation rate and litter size in sheep (Notter, 2008). Noncoding RNAs, which make up more than 98% of the genes in the human genome, can be identified and detected using RNA‐seq. Noncoding RNA transcripts with a length of more than 200 nucleotides are known as LncRNAs (Kopp & Mendell, 2018). In 1998, the first LncRNAs were identified in mice, and they were designated Human Homolog 19 (Pachnis et al., 1988). Many LncRNAs have recently been discovered in goats and sheep. Small‐tail Han sheep are a Chinese breed with hyper prolificacy but a slower growth rate. Multiple LncRNAs have been related to fecundity in Small‐tail Han sheep. (Miao et al., 2016b). In rats, the tissue‐specific expression of LncRNAs was found to be higher than that of mRNA (Wang et al., 2014). Multiple LncRNAs are knocked down in mature oocytes, resulting in larger blastocysts and higher growth rates in cows (Caballero et al., 2015). Studies into the impacts of LncRNAs on the reproductive rate have begun in various species, and there has been limited research into the impacts of LncRNAs on the reproductive rate in sheep. With the implementation of RNA‐seq, the identification of LncRNAs has become more systematic and reliable. Bioinformatics research, on the other hand, may help to improve the economic traits of farm animals by predicting the function of LncRNAs with high accuracy (especially sheep fertility) using RNA‐seq.

Shal sheep is one of the breeds of native Iranian sheep that has a high weight gain rate, and its meat is low in fat. Shal sheep, in addition to a high percentage of meat, have a high percentage of multiparous, so they can be a good option for live sheep reproduction and sheep fattening. Breeders and residents of the villages surrounding Qazvin depend on this breed for their livelihood (Hashemi & Ghavi Hossein‐Zadeh, 2020).

The size of Sangsari sheep is usually small. The tail of a Sangsari sheep is much smaller than that of other breeds of sheep, and such sheep usually have no tail or a very small tail. Due to their small size, Sangsari sheep have a great ability to walk, and due to the desert nature of their habitat, they are very resilient sheep. Multiparous in Sangsari sheep is very rare (0%), and one of the sheep of Iranian breed is considered uniparous and low fertility. This breed is one of the endangered breeds in Iran, and the Ministry of Agriculture is taking special initiatives to preserve the genetics of this breed. Improving Sangsari breed reproduction may help to conserve the breed's genetics. It will also improve profitability for breeders and be cost‐effective (Ahmadi‐hamedani et al., 2016; Miraei‐Ashtiani et al., 2007).

In this research, we surveyed LncRNAs with differential expression levels between the ovary tissue of Shal and Sangsari sheep to study the role of LncRNAs on fertility and lambing rate. After identifying and predicting the target genes, we tried to identify and predict their functions. The findings of this research will help us to further understand the function of LncRNAs in fertility. The aim of this research was to identify the LncRNAs that have differential expression between the ovaries of these two breeds, which may affect the fertility differences between the two breeds.

MATERIALS AND METHODS

Processing of animals

To carry out this research, three adult ewes of the same age (5 years old) from each breed of Shal and Sangsari were used. Shal and Sangsari ewes were selected from the farm of the Tarbiat Modares Faculty of Agriculture. The ewes were selected based on record pedigree details. According to the pedigree registered in the breeding centre of the Tarbiat Modares Faculty of Agriculture over several decades, we selected the ewes that were not related to each other. Each ewe had three consecutive lambing records. Shal ewes were multiparous per lambing, but Sangsari ewes were uniparous per lambing. All ewes were bred with unlimited access to water and food. Ewes were first ensured that they were not pregnant, and then oestrus synchronization was performed using the Gordon protocol (Gordon, 1975): Day 1, cider insertion into ewes uterus; Day 3 prostaglandin (d‐cloprostenol) injection; Day 12, cider removal and injection of PMSG; and finally day 13, slaughter and removal of ovarian tissue containing follicles being released into the fallopian tube. The ovarian tissue of these ewes was isolated and immediately transferred to the laboratory using liquid nitrogen and kept at −70°C until use for RNA extraction.

RNA extraction

A GeneJET TM RNA purification kit (Thermo Science, USA) was used to extract the total RNA according to the manufacturer's instructions. The quality and quantity of the extracted RNA were checked by electrophoresis, and RNA sequencing was carried out by BGI (China, Beijing). To avoid sequencing of ribosomal RNA, mRNA was chosen based on the poly‐A tail, and sequencing was performed using standard methods and a pair‐end library with a reading length of 150 pairs.

Bioinformatics analysis, target gene prediction and functional annotation analysis

First, FASTQC software (version 0.11.9) was used to evaluate the overall quality of the raw sequenced reads (Andrews, 2017). Then, low‐quality bases were trimmed using Trimmomatic software (version 0.4) (Bolger et al., 2014). After trimming the raw reads sequencing to ensure the post‐trimming quality improvement, quality control of the reads was performed again using FASTQC software. Reference genomes were downloaded from the Ensembl genome browser 104. Mapping of trimmed reads with the sheep reference genome (Oar_rambouillet_v1.0) was performed using STAR software (version 2.5.4a) (Dobin et al., 2013). Only reads with unique mapping were used for subsequent analyses. Samtools software (version 2.8.0) was used to sort and index the BAM file (H. Li et al., 2009). Novel LncRNAs with lengths of more than 200 nt and exon numbers greater than two were distinguished with FEELnc software after transcriptome assembly (Wucher et al., 2017).  The selection of this software was based on a previous study of our research group compering different bioinformatic tools for LncRNA predictin (Weikard et al., 2018). The FEELnc analysis consists of three phases: (1) filtering of candidate LncRNAs transcripts, (2) exploring the coding potential and nucleotide composition of candidate LncRNAs transcripts, and (3) classification of the final predicted LncRNAs transcripts. In the first phases, protein‐coding genes, monoexonic transcripts (except for those in the antisense direction to a coding gene) and transcripts with a size less than 200 NT were excluded, and a file was generated comprising potential candidate LncRNA transcripts. In the second step, the coding potential score was determined for each of the candidate LncRNA transcripts by considering the absence of an open reading frame and k‐mer nucleotide composition using the shuffle mode option of FEELnc (Chitneedi et al., 2021). The Feature Counts application (version 2.0.3), which is part of the Sub read software package, was used to create the Count matrix. LncRNAs with significant expression differences were identified in R using the DESeq2 package (version 1.32.0) (Dillies et al., 2013). Normalization of LncRNA counts was performed using the default normalization process provided by the DESeq2 package, the Median method. The threshold used to select LncRNAs with significant expression differences was and FDR < 0.05.

To investigate the possible functions of LncRNAs, we anticipated their possible target genes. The threshold was set for coding genes 100 k downstream and upstream of LncRNAs. The target genes were discovered using the biomart portion of the Ensembl genome browser 104. The target genes were analyzed using the DAVID database. The ontology terms relevant to the described differentially expressed genes were selected using a significance level (FDR < 0.05). KEGG pathways associated with genes with major expression variations were examined through the DAVID database. The KEGG pathways were selected using a significance level (FDR < 0.05). Interactions between target genes were performed using the STRING database and Cytoscape software (version 3.8.2) (Shannon et al., 2003). The degree (out degree) parameter was applied to identify hub genes using Cytoscape software.

Validation of RNA‐seq data by real‐time PCR

The expression level of several LncRNAs with their target gene was investigated to confirm the RNA‐seq results. A GeneJET TM RNA purification kit (Thermo Science) was used to extract total RNA according to the manufacturer's instructions. After determining the quality of the extracted RNA, 1 μg of RNA was used to synthesize cDNA (Step 1: for all RNA samples, 1 μg of RNA was mixed with 1 μl of Oligo (dT) primer, and then the final volume reached 13.4 μl using DEPC‐treated water. The resulting mixture was incubated for 5 min at 70°C and then immediately placed on ice. Step 2: a cDNA synthesis mixture containing 4 μl of 5× first‐strand buffer, 1 μl of dNTPs (10 mM each), 0.5 μl of RNasin (40 U/μl) and 1 μl of M‐MLV was added to each vial. Step 3: the vials were incubated for 42 min at 42°C and then for 5 min at 70°C.).

The expression level was determined by the ∆∆‐2CT method. PCR steps were performed as follows: ( (950C/5 min) 1 cycle—(950C/10 s, 580C/10 s, 720C/20 s) 35 cycles, respectively. Reverse transcription polymerase chain reaction (RT‐PCR) was repeated three independent times for each sample (Kubista et al., 2006). GAPDH was used as a reference gene. The expression levels of LncRNAs and their target genes are depicted in Figure 1. The RT‐PCR results confirmed the RNA‐seq results.

FIGURE 1

The results of the real‐time PCR. Long noncoding RNAs (LncRNAs) (ENSOARG00000026496, ENSOARG00000026495, ENSOARG00000026946, ENSOARG00000025695, ENSOARG00000025943, ENSOARG00000025639 and ENSOARG00000026864). Target genes: PLEC, CYC1, RPS9, PES1, SURF6, ITGB2 and RPL30

RESULTS

Statistical summary of sequencing and mapping the read with the reference genome

An average of 48,628,428 raw reads of 150 bp per sample were contained in the RNA‐seq data. The minimum and maximum alignments were 43,670,572 and 54,005,820, respectively. The specific alignment rate for the samples was between 78.58% and 87.26% after aligning the readings with the sheep reference genome downloaded from the Ensembl database. Table 1 provides a description of the findings obtained from mapping with the reference genome.

TABLE 1

Results of mapping with the reference genome

Sample	Shal 1	Shal 2	Shal 3	Sangsari 1	Sangsari 2	Sangsari 3
Total number of reads	48884714	59745306	48452530	44714176	44942544	45031302
Total number of alignments	45020936	54005820	43670572	47107359	46998648	46971722
Number of unique alignments	38703538	46948138	38197494	39015392	39146708	38849650
Number of readings with secondary alignment	4127902	4553838	3522662	5482386	5282062	5466856
Number of pair readings in sequencing	40893034	49451982	40147910	41624964	41716586	415048866
Percentage of unique aligned readings	79.17	78.58	78.83	87.26	87.1	86.27
Cover size	2.87	3.54	2.87	2.65	2.66	2.67

DE‐LncRNAs recognition and chromosomal classification

Only 124 of the total LncRNAs detected had different expression levels at FDR < 0.05 between Shal and Sangsari breeds (106 LncRNAs upregulated and 18 LncRNAs downregulated). The highest number of LncRNAs was found on chromosomes 4 and 17 (10 LncRNAs), while the lowest number of LncRNAs was found on chromosome 26 (without LncRNAs). The shortest LncRNAs were 731 nucleotides (ENSOARG00000025757), while the longest LncRNAs were 462337 nucleotides (ENSOARG00000026777). Figure 2 illustrates the distribution of DE‐LncRNAs on chromosomes (FDR < 0.05). Table 2 displays the top 10 LncRNAs with the highest expression levels and lowest expression levels in sheep breed Shal compared to Sangsari sheep (FDR < 0.05).

FIGURE 2

The chromosomal distribution of 124 DE‐long noncoding RNAs (LncRNAs) was significant (FDR < 0.05)

TABLE 2

Top 10 DE‐long noncoding RNAs (LncRNAs) with high and low level of expression in Shal sheep compared to Sangsari sheep (p adj < 0.05)

High expression				Low expression
Ensemble ID	Chr	log2FC 1	Padj	Ensemble ID	Chr	log2FC	p adj
ENSOARG00000026482	18	4.11	5.70E‐21	ENSOARG00000026272	10	−1.67	0.001
ENSOARG00000025779	2	4.02	1.71E‐22	ENSOARG00000025798	2	−1.40	2.62E‐09
ENSOARG00000026669	19	3.15	2.35E‐71	ENSOARG00000025540	1	−1.35	8.95E‐04
ENSOARG00000026030	3	2.13	2.30E‐15	ENSOARG00000026092	21	−1.32	0.01
ENSOARG00000025716	17	2.02	2.32E‐10	ENSOARG00000025667	17	−1.31	0.0005
ENSOARG00000026108	21	1.99	0.0002	ENSOARG00000026020	3	−1.28	2.21E‐02
ENSOARG00000026779	22	1.92	0.0002	ENSOARG00000025592	1	−1.18	3.57E‐06
ENSOARG00000025492	1	1.88	9.04E‐05	ENSOARG00000026462	18	−1.16	9.96E‐09
ENSOARG00000025946	3	1.84	2.95E‐05	ENSOARG00000025367	X	−1.16	4.43E‐11
ENSOARG00000025857	2	1.68	0.0005	ENSOARG00000025218	4	−1.13	3.86E‐02

The adjustment coefficient is shown with two decimal points. [Correction added on 21 June 2022, after first online publication: The table title was updated in this version.]

Analysis of the target genes of LncRNAs with differential expression

Ontological analysis of these LncRNAs’ target genes revealed that they were involved in 28 biological processes, 36 cellular components, 23 molecular functions, and 100 KEGG pathways. Eventually, among these, one biological process, one cellular component, and 21 KEGG pathways were significant (FDR < 0.05). The role of the target genes was calculated using the DAVID database in Figure 3.

FIGURE 3

The role of the target genes was evaluated using the DAVID database. Biological processes (red), cellular component pathways (yellow) and KEGG pathways (blue)

Interaction networks

Interaction networks of the target genes were developed using the STRING database and Cytoscape software to further expose the possible functions of target genes in the reproductive process. Figure 4 demonstrates the effects of the interaction between the target genes.

FIGURE 4

The outcomes of the interaction between the target genes as determined by the STRING database and Cytoscape software. The hub genes are depicted by green circles

DISCUSSION

The ovaries, which are one of the most important sexual organs in ewes, control follicle growth and hormone secretion (McGee & Hsueh, 2000). The capacity of an oocyte to fertilize, mature, and produce natural offspring can be characterized as its quality (Duranthon & Renard, 2001). To explore the function of LncRNAs in ovaries, we used RNA‐seq to identify the expression level of LncRNAs in the ovaries of Shal sheep and Sangsari sheep. Only 124 of the total LncRNAs detected had different expression levels at FDR < 0.05 (106 LncRNAs upregulated and 18 LncRNAs downregulated).

The PES1 gene is the target gene of ENSOARG00000025695 (LncRNA), and it has the highest degree (Figure 4). The ENSOARG00000025695 and PES1 gene have high expression levels in the Shal breed (the PES1 gene expression level was assessed by RT‐PCR) and we can assume that this correlation can result in higher fertility. The ovary is the main source of oestrogen and one of the most important target organs, since this hormone is essential for the production and function of sex tissues (Rosenfeld et al., 2001). In epithelial ovarian cancer, the PES1 gene controls the expression of oestrogen receptor alpha and oestrogen receptor beta. The loss of the PES1 gene resulted in a major downregulation of oestrogen receptor alpha and oestrogen‐responsive genes, including HIF1α, cyclin D1 and VEGF, as well as an increase in p21WAF1 and oestrogen receptor beta. In line with our findings, PES1 gene expression was found to be positively associated with oestrogen receptor alpha expression and negatively correlated with oestrogen receptor beta expression (W. Li et al., 2013).

The ribosomal protein S9 (RPS9) gene is located on chromosome 14 and is a target gene of ENSOARG00000026946 (LncRNAs). The RPS9 gene participated in the ribosome pathway (KEGG pathway: oas03010). The RPS9 gene and ENSOARG00000026946 have high expression levels in the Shal breed (RPS9 gene expression level was assessed by RT‐PCR), and we can assume that this correlation can result in higher fertility, as in the Australian dairy cow population, RPS9 gene variants (one upstream, one intron, and six downstream) had the greatest association with fertility (S. G. Moore et al., 2016).

The EF1 gene was the target gene of ENSOARG00000026495 and ENSOARG00000026496. The EF‐1 gene participated in the herpes simplex infection pathway (KEGG pathway: oas05168). The EF‐1 gene had a high level of expression in the Shal breed (EF‐1 gene expression level was assessed by RT‐PCR). The EF1 gene has the highest expression in oocytes and blastocysts, according to an analysis of 536 housekeeping genes in humans (Kakourou et al., 2013). Plectin is one of the most abundant polypeptides in Chinese hamster ovary cells, accounting for up to 1% of total protein. It is also one of the most effective phosphorus receptors (Herrmann & Wiche, 1983).

The Plectin genes were the target genes of ENSOARG00000026495 and ENSOARG00000026496. ENSOARG00000026496 had a common target gene of Plectin and EF‐1; on the other hand, these two genes have a significant and high expression level in Shal sheep (Plectin gene expression level was assessed by RT–PCR). ENSOARG00000026495 and ENSOARG00000026496 have high levels of expression in the Shal breed. Figure 4 indicates that the RPS9 gene has been shown to have a direct connection with both the PeS1 and Plectin genes. The EF‐1 gene, on the other hand, has a direct connection with the Plectin gene. Since these genes have a straight connection, they are likely to function together to influence fertility. Finally, we can assume that the correlation of these genes together with ENSOARG00000026495 and ENSOARG00000026496 can result in higher fertility. The PPP1R42 gene is an ENSOARG00000026531 target gene. The PPP1R42 gene participated in molecular function (GO: 0005813). Shal sheep had higher expression of ENSOARG00000026531 than Sangsari sheep. The PPP1R42 gene was one of 50 endometrial genes (upregulated) found to be a biological marker for determining the pregnancy status in cows on the seventh day of the oestrous cycle (Rabaglino & Kadarmideen, 2020).

The EEF1A2 gene is an ENSOARG00000026240 target gene. The EEF1A2 gene participated in GTPase activity (GO: 0003924) and GTP binding (GO: 0005525). In the late embryonic phase in sheep, the EEF1A2 gene is positively associated with embryo weight and plays an important function in muscle growth during the embryonic phase (Xu et al., 2014).

The MAPK1 gene is a target gene of ENSOARG00000025710 and ENSOARG00000025711. The MAPK1 gene participated in the oestrogen signalling pathway (KEGG pathway: oas04915), cAMP signalling pathway (KEGG pathway: oas04024) and MAPK signalling pathway (KEGG pathway: oas04010). The MAPK1 gene is essential for the maturation and growth of oocytes (López‐Cardona et al., 2017). MAPK gene expression levels over a certain threshold are a significant predictor of ovule puberty and a criterion for ovule puberty efficiency (Y. Li et al., 2017; Sun et al., 2016). MAPK gene active molecules have also been shown to positively pass BMP signals (R. K. Moore et al., 2003). In mice, MAPK1 gene activity is still critical for the performance of suitable granulosa cells prior to ovule puberty (Gratao et al., 2008).

PRKACA (protein kinase catalytic subunit alpha) is a target gene of ENSOARG00000025288. The PRKACA gene participated in the oestrogen signalling pathway (KEGG pathway: oas04915), cAMP signalling pathway (KEGG pathway: oas04024) and MAPK signalling pathway (KEGG pathway: oas04010). The ENSOARG00000025288 and PRKACA gene have high expression levels in the Shal breed, and we can assume that this correlation can result in higher fertility. The insulin signalling pathway's network partners PRKACA, CALM and PPP1CC have been identified as controlling the resumption of oocyte meiosis. Several regulatory proteins are dephosphorylated by PPP1CC during cell division. The active form of PPP1CC dephosphorylates calcium/calmodulin‐dependent protein kinases during insulin signalling, promoting glucose metabolism. The activation of PRKACA is activated by the interaction of INS or IGF1 with its recipient on the oocyte (Selvaraju et al., 2018). The PRKACA gene plays a role in fertility‐related signalling pathways such as oxytocin, gonadotropin‐releasing hormone, oestrogen, and transforming growth factor beta (La et al., 2019). The function of mitochondria in oocyte maturity, fertilization, and embryonic growth has long been understood (Van Blerkom et al., 1995). Since an oocyte must have an optimum number of viable mitochondria to be of high growing ability, reducing the number of viable mitochondria in the oocyte can have a negative effect on embryonic growth. The ability to generate ATP is essential for nuclear and cytoplasmic puberty in preparation for fertilization and completion of meiosis II with adequate mitochondrial numbers and ATP levels (John et al., 2010). Reduced ATP production and the number of mtDNA copies are linked to poor oocyte quality and reduced fetus growth (Wakefield et al., 2011). The hub genes in Figure 4 included mitochondrial genes such as CYC1, MRPL12, MTG1 and MRM2. The cytochrome c1 gene (CYC1) is the target gene of ENSOARG00000026495 and ENSOARG00000026496. The CYC1 gene participated in metabolic pathways (KEGG pathway: oas01100). The CYC1 gene with Plectin was a common target gene of ENSOARG00000026495 and ENSOARG00000026496, and CYC1 with EF1 genes was a common target gene of ENSOARG00000026496. In the Shal breed, the CYC1, Plectin and EF1 genes have high expression levels (the CYC1, Plectin and EF1 gene expression levels were assessed by RT‐PCR), and we can conclude that this association leads to increased fertility. As a result, ENSOARG00000026496 can influence fertility through the CYC1 and EF1 genes, and ENSOARG00000026495 and ENSOARG00000026496 can influence fertility through the CYC1 and Plectin genes. The CYC1 gene is an essential component of an apoptotic pathway that reacts to DNA damage and other forms of cellular stress. In fact, mouse fetuses missing the CYC1 gene die in the womb in the middle of pregnancy (K. Li et al., 2000). On the other hand, uteruses that developed a living fetus have a high level of CYC1 protein, indicating that this protein is essential for fetus survival (Beltman et al., 2014).

The BRD2 gene was found to be directly related to the PES1 gene, as shown in Figure 4. The BRD2 gene is a target gene of ENSOARG00000026854. The BRD2 gene participated in chromatin binding (GO: 0003682). The BRD2 gene is one of the genes with extreme expression in the testes and ovaries (Rhee et al., 1998). BRD2 interacts with the GDF9 and BMP15 genes in the protein–protein interaction network of female mammals (Ahmad et al., 2017).

The integrin subunit beta 2 (ITGB2) gene is the target gene of ENSOARG00000025639. The ITGB2 gene was found to be directly related to the Plectin gene (Figure 4). Recent research proposes that integrin could carry out a function in Gamete‐special interactions and mergers throughout fertilization in pigs (Linfor & Berger, 2000). Oocyte morphology and consistency improved as ITGB2 protein levels increased (Antosik et al., 2010), but in our research, there was no significant difference in ITGB2 gene expression levels between the two breeds (the ITGB2 gene expression level was assessed by RT‐PCR).

The SURF6 gene was the target gene of ENSOARG00000025943. In the Shal breed, the SURF6 gene and ENSOARG00000025943 have high expression levels, and we can conclude that this association leads to increased fertility. The SURF6 gene is located in the nucleus of transcription‐capable cells and has a specific distribution in mouse oocytes and preimplanted fetuses (Romanova et al., 2006).

CONCLUSIONS

This is the first study to use RNA‐seq data to classify LncRNAs, GO terms and interactions of LncRNA target genes in ovarian tissue of Iranian sheep. We discovered an important association between the target genes of the LncRNAs identified in our research. The LncRNAs that are uniquely enabled in sheep prolificacy are systematically revealed in this analysis. Large‐scale LncRNAs studies exploring the molecular mechanisms of sheep fecundity and finding the related LncRNAs leading to fecundity variations will be made simpler with this approach. Our results include a strategy for identifying possible candidate LncRNAs as well as a blueprint for future efforts to address sheep prolificacy barriers.

CONFLICT OF INTEREST

The authors declare no conflict of interest.

AUTHOR CONTRIBUTIONS

Data curation, formal analysis, methodology, software, validation, writing—original draft and writing—review and editing: Shahram Hosseinzadeh. Data curation, formal analysis, methodology, project administration, supervision and writing—original draft: Ali Akbar Masoudi. Formal analysis and software: Alireza Ehsani and Rasoul Vaez Torshizi.

ETHICS STATEMENT

All procedures involving animals were approved by the scientific committee of the Department of Animal Science, Tarbiat Modares University.

PEER REVIEW

The peer review history for this article is available at https://publons.com/publon/10.1002/vms3.859.