Identification of non-coding and coding RNAs in porcine endometrium

    Identification of non-coding and coding RNAs in porcine endometrium Yueying Wang, Tao Hu, Lihang Wu, Xiaoran Liu, Songyi Xue, Minggang Lei PII: DOI: Reference:

S0888-7543(16)30127-6 doi: 10.1016/j.ygeno.2016.11.007 YGENO 8850

To appear in:

Genomics

Received date: Revised date: Accepted date:

25 June 2016 16 November 2016 26 November 2016

Please cite this article as: Yueying Wang, Tao Hu, Lihang Wu, Xiaoran Liu, Songyi Xue, Minggang Lei, Identification of non-coding and coding RNAs in porcine endometrium, Genomics (2016), doi: 10.1016/j.ygeno.2016.11.007

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

Identification of non-coding and coding RNAs in porcine

T

endometrium

1

SC R

Lei1*

IP

Yueying Wang1 · Tao Hu1 · Lihang Wu1 · Xiaoran Liu1 · Songyi Xue1 · Minggang

Key Lab of Swine Genetics and Breeding of Ministry of Agriculture, Key Lab of

NU

Agricultural Animal Genetics, Breeding and Reproduction of Ministry of Education, College of Animal Science and Technology, Huazhong Agricultural University,

MA

Wuhan, Hubei, PR China *

D

Correspondence: Minggang Lei, College of Animal Science and Technology,

TE

Huazhong Agricultural University, Wuhan, 430070, Hubei, PR China. Email:

AC

CE P

[email protected]

Page 1

ACCEPTED MANUSCRIPT Abstract One of the most critical periods of embryonic loss in pig is Day 12 of pregnancy, when implantation begins. Here, we analyzed the gene expression on Day 12 of pregnancy and non-pregnancy in the porcine endometrium using RNA

IP

T

sequencing (RNA-seq). 237 mRNAs, 34 lncRNAs and 1 miRNA were significantly differentially expressed between the two groups. Further functional analyses were

SC R

conducted to identify these differentially expressed transcripts. The results demonstrated that they participate in various biological processes, such as cell adhesion, binding, nucleic and metabolic processes. In addition, our results showed

NU

that the differentially expressed genes (IL1R, FGF9, DUPS10, DUPS4, CD14 and MAP4K4) in MAPK pathway, and lncRNAs of XLOC_2604764 and XLOC_2604756

MA

may play a vital role in regulating embryo implantation. Besides, we investigated the lncRNA-ssc-miR-132-mRNA interactions, aiming to explain the regulatory networks

D

of coding and non-coding genes that contributes to the establishment of the maternal

Keywords

Embryo

TE

pregnancy.

implantation

·

Non-coding

RNA

·

Coding

AC

Introduction

CE P

RNA· Pig · Endometrium

One of the most critical periods of embryonic loss is Day 12 of pregnancy, when the maternal recognition of pregnancy and implantation begin [1, 2]. Successful implantation requires the activity of conceptus and dramatic morphological and functional changes in the uterine endometrium to make the uterus receptive to the developing conceptus; besides, the conceptus produces estrogen, which acts as a signal for maternal recognition of pregnancy in pigs [3, 4]. Accordingly, understanding the mechanisms of endometrium changes during implantation process is essential to increasing the rate of successful pregnancies. Implantation is a complex and well-coordinated interaction between the developing Page 2

ACCEPTED MANUSCRIPT conceptus and the maternal uterus [5]. On Day 12 of pregnancy, embryonic growth and development are dependent on hormones, growth factors and transport proteins in the uterine lumen, all of which are secreted primarily by the epithelia of maternal

IP

T

uterus [6-8]. In turn, the uterine endometrium becomes receptive to the conceptus through the regulation by cell adhesion molecules, growth factors and cytokines [9,

SC R

10]. Previous studies have reported that the genes in endometrium have different expression patterns between pregnant and non-pregnant pigs, such as the growth factors (FGF7 [11, 12], EGF [13, 14], HB-EGF [15, 16] and TGF-α [17]), adhesion

NU

molecules (L-selectin [18], IFN [4, 19] and osteopontin [20, 21]), and cytokines (LIF [22, 23], IL-6 family [24] and IL-11 [25]). Besides, some miRNAs have also been

MA

reported to regulate embryo implantation process via regulating the expression of particular genes. Recent study has revealed distinctive miRNAs expression patterns in

D

porcine endometrium during implantation, and suggested that miR-181a and

TE

miR-181c likely play important roles in regulating the genes and pathways that are known to be involved in embryo implantation and placentation in pigs [26].

CE P

Cyclooxygenase-2, a gene critical for implantation, is post-transcriptionally regulated by two miRNAs, namely mmu-miR-199a and mmu-miR-101a [27]. All these studies were merely based on one single transcript level, namely mRNA or miRNA.

AC

Recently, a number of long non-coding RNAs (lncRNAs) have been shown to compete with mRNAs for binding to miRNAs and contribute to many biological processes. Non-coding RNAs (ncRNAs) are roughly grouped into two important types based on their transcript sizes: small ncRNAs and lncRNAs. MiRNAs are small RNA molecules that bind to the 3’-untranslated region (3’-UTR) of mRNAs, and regulate gene expression via translational inhibition or degradation of mRNA transcripts. The transcripts of larger ncRNA molecules longer than 200 bp can be classified into lncRNAs, which have been shown to play a role in multiple cellular maintenance functions, such as protein scaffolding, chromatin looping, and regulation of mRNA stability [28, 29]. There are growing evidences showing that lncRNAs may Page 3

ACCEPTED MANUSCRIPT potentially interact with miRNAs and modulate each other’s expression: (i) lncRNAs may act as competing endogenous RNAs (ceRNAs), namely miRNA sponges or antagonist, which may down-regulate the expression and activities of miRNAs, and

IP

T

subsequently regulate the expression of miRNA target genes at post-transcriptional level; and (ii) miRNAs can inhibit the expression of lncRNAs [30-32]. For example,

SC R

as for human urothelial carcinoma associated 1 (UCA1), a previous study has shown that lncRNA UCA1 may function as an endogenous sponge to up-regulate the expression of fibroblast growth factor 1 (FGFR1) through directly binding and

NU

inhibiting the expression of miR-216b, and UCA1 is also involved in the activation of ERK signaling pathway [33].

MA

In the present study, we analyzed the gene expression in porcine endometrium on Day 12 of pregnancy and non-pregnancy using next-generation RNA sequencing

D

(RNA-seq). To our knowledge, no RNA-seq data that demonstrate the molecular

TE

changes in Duroc pig endometrium are available so far. This study was aimed to obtain deep insights into endometrial coding and non-coding RNA changes associated

CE P

with conceptus implantation and to identify the specific mechanisms involved in the regulation of implantation. Furthermore, one advantage of this work is that we investigated the lncRNA-mRNA-miRNA cross-talk that mediates the conceptus

AC

implantation in endometrium of pig for the first time. The findings of the present study are expected to facilitate the understanding of the mechanisms that underlie successful conceptus implantation in pig. Materials and methods Sample collection Six Duroc gilts with similar age and genetic background from one commercial herd were selected. Animals were randomly assigned to non-pregnancy (n=3) (no mating estrus) and pregnant (n=3) groups. Gilts of pregnant group were artificially inseminated twice after estrus. Uteri were obtained from animals slaughtered on Day Page 4

ACCEPTED MANUSCRIPT 12 (n=3) of both groups. Each uterine horn was flushed with sterile PBS (pH 7.4), and subsequently opened longitudinally at the inner side. The formation of the conceptus on day 12 of pregnancy was filamentous. Samples from the endometrium of the

IP

SC R

nitrogen and stored at -80°C before RNA isolation.

T

pregnant and non-pregnant sows were taken. Tissue samples were frozen in liquid

Total RNA isolation

Total RNA was isolated from each individual sample using TRIzol reagent (Invitrogen

NU

Corp., CA, USA). Purity and quantity of total RNA were measured by using Nanodrop equipment (Thermo Fisher Scientific, Wilmington, USA). Integrity of RNA

MA

was assessed using the RNA Nano6000 Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies, CA, USA).

TE

D

Library preparation for lncRNA sequencing and data analysis A total of 3 μg RNA per sample was used as input material for the RNA sample

CE P

preparations. Firstly, ribosomal RNA was removed by Epicentre Ribo-zero™ rRNA Removal Kit (Epicentre, USA), and rRNA free residue was cleaned up by ethanol precipitation. Subsequently, sequencing libraries were generated using the

AC

rRNA-depleted RNA by NEBNext® Ultra™ Directional RNA Library Prep Kit for Illumina® (NEB, USA) following the manufacturer’s recommendations. In order to select cDNA fragments of the preferred 150~200 bp in length, the library fragments were purified with AMPure XP system (Beckman Coulter, Beverly, USA). The libraries were sequenced on an Illumina Hiseq 2500 platform and 100 bp paired-end reads were generated. Raw reads of fastq format were firstly processed through in-house perl scripts. In this step, clean data (clean reads) were obtained by removing the reads that contained adapter or ploy-N and low-quality reads from the raw data. The obtained reads were mapped with Tophat (v2.0.9) to the porcine genome sequence assembly (Sscrofa 10.2). The mapped reads of each sample were assembled Page 5

ACCEPTED MANUSCRIPT by both Scripture (beta2) [34] and Cufflinks (v2.1.1) [35]. CNCI (v2) [36] (Coding-Non-Coding-Index), CPC (0.9-r2) [37] (Coding Potential Calculator), Pfam-scan (v1.3) [38] and phyloCSF (v20121028) [39] (phylogenetic codon

IP

T

substitution frequency) were used to distinguish mRNA from lncRNAs (Fig 1). The phast software (v1.3) was generally used for conservation analysis [40].

SC R

Transcripts without coding potential were our candidate set of lncRNAs. Then, we searched coding genes 10 k/100 k upstream and downstream of lncRNA as the cis target gene. Trans role of lncRNA is to be identified by the correlation of their

NU

expression levels.

MA

Library preparation for microRNA sequencing and data analysis A total of 3 μg total RNA per sample was used as input material for the small RNA

D

library. Sequencing libraries were generated using NEBNext® Multiplex Small RNA

TE

Library Prep Set for Illumina® (NEB, USA.) and index codes were added to assign sequences to each sample. Briefly, RNAs were ligated to 3' and 5' adaptors

CE P

sequentially, reverse transcribed to cDNA and then PCR amplified. PCR products were purified on a 8% polyacrylamide gel (100V, 80 min). DNA fragments corresponding to 140~160 bp (the length of small noncoding RNA plus the 3' and 5'

AC

adaptors) were recovered and dissolved in 8 μL elution buffer. The clustering of the index-coded samples was performed on a cBot Cluster Generation System using TruSeq SR Cluster Kit v3-cBot-HS (Illumia) according to the manufacturer’s instructions. After cluster generation, the library preparations were sequenced on an Illumina Hiseq 2500 platform (Illumina, San Diego, CA, USA). Raw data (raw reads) of fastq format were firstly processed through custom perl and python scripts. In this step, clean data (clean reads) were obtained by removing the redundant reads. The small

RNA tags

of

clean

reads

were

mapped

to

reference

sequence

(ftp://ftp.ncbi.nlm.nih.gov /genomes/Sus_scrofa) by Bowtie [41]. Modified software mideep2 and srna-tools-cli were used to search for known miRNA by matching with Page 6

ACCEPTED MANUSCRIPT miRBase21 [42]. Prediction of target gene of miRNA was performed by miRanda [43]. The numbers of miRNA reads were normalised by transcripts per million (TPM) values (TPM= (miRNA total reads/total clean reads) ×106). The DEmiRNAs between

IP

T

samples were identified by DESeq. The P-values was adjusted using the Benjamini & Hochberg method. Corrected P-value of 0.05 was set as the threshold for significantly

SC R

differential expression by default. GO and KEGG enrichment analysis

NU

GO term and KEGG pathway analyses of coding genes were performed using DAVID bioinformatics resources. We used GOSeq and KOBAS softwares to test the statistical

MA

enrichment of the target gene candidates of differentially expressed lncRNAs and miRNAs in GO terms and KEGG pathways, respectively.

TE

D

Construction of Regulatory Network

For construction of lncRNA-miRNA network, we first used Blast to remove the

CE P

precursor of miRNAs which were highly matched with lncRNAs, then used miRanda software to predict the miRNA-targeted lncRNAs (alignment score > 160 and minimum free energy of -20 kcal/mol). The coding and non-coding RNA interactions

AC

were constructed using Cytoscape v3.2.1[44]. Quantitative PCR

RNA samples from the 6 animals used for the RNA-seq experiment were analyzed by qPCR. Total cDNA was synthesized using reverse transcriptase Kit (TaKaRa, Dalian, China). qPCR was performed using LightCycler 480Ⅱ Real-Time PCR System and SYBR® Green PCR Master Mix (TaKaRa, Dalian, China). Each real-time RT-PCR reaction (in 25 μL) involved 12.5μL 2×SYBR Green Real-time PCR Master Mix (TaKaRa, Dalian, China), 1 μL of each primer, 2 μL cDNA and 8.5μL H2O. The cycling conditions included an initial single cycle (95°C for 1 min), followed by 40 Page 7

ACCEPTED MANUSCRIPT cycles of 95°C for 15 s, 60°C for 15 s, 72°C for 20 s. All amplifications were followed by dissociation curve analysis of the amplified products. Specific primers were designed using the NCBI, specificities were confirmed with BLAST and gene

(−ΔΔCt)

IP

using 2

T

expression levels were normalized with RPS20 to attain the relative expression by value methods. The statistical difference in gene expression was

SC R

analyzed by SAS during different endometrium development stages of pregnancy. The correlation between the results of RNA-seq and qPCR was calculated using correlation test.

NU

Results

MA

RNA-Seq analysis of coding RNA and non-coding RNA (lncRNA and miRNA) in endometrial tissue samples

D

We compared the gene expression in the endometrium of pregnant Duroc sows with

TE

that of non-pregnant ones. The number of reads and transcripts were shown in Table 1. Finally, 203 differentially expressed mRNA (DEmRNA) were obtained, among which

CE P

170 had higher expression and 33 had lower expression in the samples on Day 12 of pregnancy (D12) than those on Day 12 of non-pregnancy (D12K) (Table S1). 1507 novel transcripts were obtained (Table S2). 1,490 lncRNAs (mainly including

AC

lincRNA, inronicRNA and anti-senseRNA) were retained after the basic filtering and removing the potential coding transcripts (Fig. 2, Table S3). The results showed that 34 lncRNAs were significantly differentially expressed between the two groups, among which 29 genes were up-regulated and 5 were down-regulated in D12 compared with in D12K (Table S1). In addition, FGG (fibrinogen gamma chain) was the predicted target gene of two significantly differentially expressed lncRNAs (DElncRNAs),

namely XLOC_2604764

(6-fold

up-regulated

in

D12) and

XLOC_2604756 (8-fold up-regulated in D12). These results were validated by qPCR (Table S4). Notably, only one differentially expressed miRNA (DEmiRNA, ssc-miR-132) was found (Table S1). Subsequently, 10 DEGs (7 mRNAs and 3 Page 8

ACCEPTED MANUSCRIPT lncRNAs) and 1 miRNAs were randomly selected to validate the reliability of the RNA-seq results. For these genes, quantitative PCR (qPCR) was performed to

T

confirm the results obtained with RNA-seq (Table 2).

IP

Relationship between lncRNAs and their neighboring coding genes

SC R

Previous studies have reported that some lncRNAs may act in cis and affect the gene expression of their chromosomal neighbourhood [45, 46], and most lncRNA transcripts can also be derived from divergent transcription [47, 48]. To validate this

NU

hypothesis, firstly, we annotated the location relationship between each lncRNA and its cis target genes (Table S5). Secondly, we analysed the gene pairs formed by

MA

lncRNAs and their neighboring genes and identified 881 lncRNAs: coding gene pairs in total (164 in divergence) (Fig. 3a). Then, functional analysis was performed on

D

these coding genes near the predicted lncRNAs. The results indicated that “ion

TE

binding” and “cation binding” were the enriched Gene Ontology (GO) terms (Fig. S1), and “Ubiquitin mediated proteolysis”, “Insulin signalling pathway” and “Leukocyte

CE P

transendothelial migration” were the enriched pathways (Fig. S2), all of which have been confirmed to participate in the process of embryo implantation. Thirdly, the correlation between the lncRNAs and their neighbouring protein-coding genes was

AC

higher than that between the neighbouring protein-coding genes (Fig. 3b, c). Indeed, the expression patterns of lncRNAs were more strongly correlated with those of their neighbouring coding RNA (mean correlation: 0.32) compared with the correlation between the expression patterns of coding gene pairs (mean correlation: 0.26) (P-value= 0.0007152, Kolomogorv [KS] test) in all neighbouring pairs. However, lncRNA: coding gene pairs exhibited similar expression patterns (mean correlation: 0.26) with coding gene pairs (mean correlation: 0.31) in divergent, even though there was a slight difference in their mean correlations (P-value= 0.1579, [KS] test). Furthermore, the correlation of lncRNA: coding gene pairs was higher than that of random ones in all pairs (p-value= 2.552e-06, [KS] test), especially in divergent pairs Page 9

ACCEPTED MANUSCRIPT (p-value < 2.2e-16, [KS] test). Fourthly, it was shown that the 5’ ends of lncRNAs were enriched in a 4 kb region surrounding the transcription start sites (TSSs) of their neighbouring coding genes (Fig. 3d), which is consistent with the results of previous

IP

T

studies [49]. Previous studies have shown that divergent transcription is an important source of lncRNAs [48, 50]. These results also reveal that most lncRNA are

SC R

originated from TSS regions of the coding genes in a divergent way. Functional analysis of differentially expressed genes in endometrial tissue samples

NU

To evaluate the biological functions and signal transduction pathways regulated by pregnancy in the uterine endometrium, GO and KEGG analysis were performed on

MA

the obtained mRNAs, lncRNAs and miRNAs. In accordance with previous studies [51, 52], GO terms of the up-regulated DEmRNAs were significantly enriched in the

D

biological processes of mitosis, epidermis development and cell adhesion and other

TE

four molecular function terms throughout DAVID program (Table 3). “ABC transporters”, “MAPK signalling pathway” and “regulation of actin cytoskeleton”

CE P

were the only three enriched pathways (Table 3). Also, “MAPK signalling pathway” is a significant pathway for the successful implantation of embryo. DGEs IL1R, FGF9, DUPS10, DUPS4, CD14 and MAP4K4, which may be involved in the regulation of

AC

embryo implantation, were in the pathway of MAPK. To further predict the regulatory function of lncRNAs in endometrium of pigs, we performed GO analysis for the target genes predicted with DElncRNAs in cis and trans by GOSeq. As the number of DElncRNAs was small, there were no GO terms significantly enriched for cis target genes, and some GO terms for trans target genes were enriched in 16 processes, such as “nucleic and metabolic process”, “cellular nitrogen compound metabolic process” and “binding” (Fig. 4a). According to the KEGG database, the most significantly enriched pathway was “cell cycle”. Regarding the DEmiRNA (ssc-miR-132), 10 GO terms were significantly enriched for the target genes, such as “cation binding” and “transition metal ion binding” (Fig. 4b). Page 10

ACCEPTED MANUSCRIPT KEGG analysis showed that the most enriched pathways were “cell cycle” and “adherens junction”. Interestingly, the most significantly enriched GO terms and pathways were similar between the target genes of DElncRNAs and DEmiRNA,

IP

T

indicating that they may have similar regulation functions with mRNAs or may

SC R

regulate each other at the phase of implantation. Construction of regulatory network

Interactions of lncRNA-miRNA were constructed. As miR-132 was the only

NU

DEmiRNA, we investigated the interaction of miR-132 with the predicted lncRNAs and mRNAs. The miRanda software identified 6 lncRNAs (XLOC_2558587, XLOC_1796741,

XLOC_2087578,

MA

XLOC_782816,

XLOC_933841

and

XLOC_532476) with potential binding sites of miR-132. Therefore, TIGD2 (tigger

D

transposable element derived 2) was the predicted target gene of a lncRNA in cis,

TE

namely XLOC_2558587. GDF2 (growth differentiation factor 2), B3GAT1 (beta-1, 3-glucuronyltransferase 1), LOC100515931, LOC100614836 and IL15RA (interleukin

CE P

15 receptor, alpha) were the predicted target genes of lncRNA XLOC_782816, XLOC_1796741, XLOC_2087578, XLOC_933841 and XLOC_532476 in trans,

AC

respectively. Discussion

Establishment of the cross-talk between the developing conceptus and the maternal endometrium involves a complex process, particularly when the conceptus has already been elongated and begins to be attached to the uterine endometrium around day 12 of gestation [53]. At this time, gene expression in endometrium will be changed for the attachment of conceptus to the endometrial epithelial cells [54]. To better understanding the implantation process, we profiled the endometrial coding and non-coding RNA expression and identified the genes differentially expressed in D12 and D12K endometrium by using RNA-seq. Numerous studies have been conducted Page 11

ACCEPTED MANUSCRIPT to detect the molecular mechanisms of implantation in pigs. Zhang et al. (2013) profiled the gene expression in endometrium on day 12 of gestation and selected 13612 DEGs between the two breeds using digital gene expression profiling [55].

IP

T

Microarrays have been used to describe endometrial gene expression on day 12 of the estrous cycle and pregnancy [56-58]. Similarly, systematic studies have been

SC R

performed by RNA-seq to describe the transcriptome changes in prepubertal gilt (crossbreed of German Landrace and Pitrain) endometrium on day 12 and 14 during estrous cycle and pregnancy, and 2593 DEGs were revealed from the comparison of

NU

day 12 pregnant animals and corresponding non-pregnant control [51, 52].

And the

main difference is that all of these studies focus on the coding RNAs in Yorkshire,

MA

Meishan or crossbreed pigs. In this study, we analyzed the coding and non-coding RNA expression in Duroc porcine endometrium. As for MAPK pathway, we enriched

D

the pathway in this study which was not mentioned in the previous study.

TE

In this study, a number of genes and transcripts were detected to be differentially expressed in endometrium tissue between the two groups. Bioinformatics analysis

CE P

showed that the up-regulated coding genes are mainly involved in cell-to-cell interaction and tissue remodeling, and are enriched in MAPK pathway genes. Furthermore, functional clustering analysis of the target genes of DElncRNAs and

AC

DEmiRNAs showed that “cell cycle” and “adherens junction” were the shared GO terms. 6 DGEs, namely IL1R2, FGF9, DUPS10, DUPS4, CD14 and MAP4K4, are involved in the MAPK signalling pathway. Some of these genes contribute to the establishment of maternal-fetal immune tolerance, such as IL1R2 and CD14. Interleukins play a major role in many biological processes including inflammation response or adipogenesis, but they also seem to be implicated in critical reproductive events such as proliferation and differentiation of cells and participate in embryo-maternal dialog [59]. Additionally, IL1R2 is the receptor of IL-1, and can capture IL-1 and inhibit its binding to another receptor IL-1R1, suggesting that IL1R2 plays an important role in regulating the biological activities of IL-1 [60, 61]. The Page 12

ACCEPTED MANUSCRIPT reduced expression of IL1R2 within the implantation window suggests the existence of accurate regulatory mechanisms that alleviate IL-1 inhibition during this crucial period and facilitate IL-1 pre-implantation actions by down-regulating IL1R2

IP

T

expression in human [62]. In contrast, the expression of IL1R2 in D12 was 8-fold higher than in D12K, which might be attributed to the different implantation

SC R

processes between human and pig.

Many genes, which regulated the embryo implantation, were up-regulated in D12, such as FGF9, IGFBP2 and ESR2 (3-, 4-, 5-fold up-regulated in D12, respectively)

NU

has been shown to be regulated by E2 and progesterone (Table S6) [63, 64]. Previous studies have reported that FGF9 is significantly highly expressed in pregnant animals

MA

and was identified as a growth factor in human endometrium [63]. In pig,

the

strongest staining of FGF9 was observed at the apical domain of the glandular

D

epithelial cells, but the staining of its receptor FGFR3 was significantly lower in the

TE

pregnant sows [63, 65]. Hence, it was speculated that FGF9 might not be related to the physiological and morphological changes within the endometrium, but rather be

CE P

implicated in the embryo-maternal communication [65]. Up-regulation of IGFBP2 was also found in porcine endometrium at Day 14 and 18 of pregnancy in bovine endometrium [65, 66]. Moreover, a number of DEGs were implicated in the

AC

establishment of pregnancy in pigs through calcium ion binding and regulation of calcium homeostasis, such as S100A8, S100A12 and TRPV6 (11-, 9-, 4-fold up-regulated in D12, respectively). S100G is known to be involved in the regulation of intracellular calcium levels [67], and S100A8 and S100A12 are calcium-binding proteins that belong to the S100 protein family [68]. A previous study has demonstrated that LPS causes a rapid up-regulation of S100A8, S100A9 and S100A12 released from damaged endometrial cells and may play a major antimicrobial role [69]. Ka et al. (2009) found that the calcium-related genes S100A7A, GSN and TRPV6 were expressed at higher levels in pregnant uterine endometrium [70]. TRPV6, a member of the TRPV family, acts as an epithelial calcium ion channel, Page 13

ACCEPTED MANUSCRIPT and its expression was predominantly localized to endometrial luminal epithelium (LE) cells and up-regulated in the endometrium of pregnant pig uterus [67, 70]. Hence, LE cells utilize TRPV6 to absorb calcium ions present in the uterine lumen during

IP

T

implantation period [1, 71]. But the detailed mechanisms by which S100A8 and S100A12 regulate implantation should be further analyzed.

SC R

Successful implantation in pig requires the participation of connexin proteins such as MMPs, which are endopeptidases that play an important role in the transient invasive property of trophoblastic cells during implantation[72]. Our results indicated that

NU

among these MMPs, MMP7 was up-regulated (by 7-fold in D12) and involved in Wnt signal pathway. Expression of MMP7 was extremely high in the epithelial areas as

MA

confirmed in secretory human endometrium by using microarrays [73]. Another study reported that MMP7 was expressed during the receptive phase, although to a small

D

extent, thus it is necessary to determine the role of MMP7 in endometrial receptivity

TE

[74].

Recent reports have suggested that lncRNAs can potentially interact with other

CE P

classes of non-coding RNAs including miRNAs and then regulate the expression of mRNA. In this study, we systematically analyzed the complex effects of the interactions

between

miRNAs

and

their

target

genes

and

provided

AC

lncRNA-miRNA-mRNA networks. Only 1 miRNA (miR-132) was up-regulated in D12 as compared with D12K, suggesting that the regulatory effect of miRNAs on implantation is weak. MiR-132 is involved in the cAMP signaling pathway and promotes estradiol synthesis in ovarian granulosa cells [75]. It was also observed that overexpression of miR-132 increased E2 secretion from KGN, a steroidogenic human granulosa-like tumor cell line [76]. Therefore, we speculated that miR-132 might participate in the process of implantation. The lncRNA-miR-132-mRNA networks revealed a set of lncRNAs and mRNAs that are possibly involved in embryo implantation process. This study provides new insights into the molecular mechanism of porcine embryo implantation. Page 14

ACCEPTED MANUSCRIPT Conclusions In this study, we profiled the expression of lncRNAs, miRNAs and mRNAs in pig

T

endometrium on Day 12 of pregnancy and non-pregnancy in Duroc pigs.

IP

Bioinformatics analyses suggest that some lncRNAs are involved in important

SC R

biological processes associated with embryo implantation. Moreover, we present the transcriptome-scale study on coding and non-coding RNA interactions in porcine

NU

endometrium. This study also provides valuable information for lncRNA studies in

MA

other animals with noninvasive implantation. Compliance with ethical standards

Funding This work was supported by the National Basic Research Program of China

D

(2012CB138504) and the National Porcine Industry Technology System (CAR-36).

TE

Conflict of interest The authors declare that they have no competing interests. Ethical approval All studies involving animals were conducted according to the

CE P

regulation approved by the Standing Committee of Hubei People’s Congress, P. R. China. Sample collection was approved by the ethics committee of Huazhong Agricultural University. Animals were humanely sacrificed as necessary to ameliorate

AC

suffering.

References [1] Geisert RD, Renegar RH, Thatcher WW, Roberts RM, Bazer FW. Establishment of pregnancy in the pig: I. Interrelationships between preimplantation development of the pig blastocyst and uterine endometrial secretions, Biol Reprod. 27 (1982) 925-39. [2] Bazer FW, Thatcher WW. Theory of maternal recognition of pregnancy in swine based on estrogen controlled endocrine versus exocrine secretion of prostaglandin F2alpha by the uterine endometrium, Prostaglandins.14 (1977) 397-400. [3] Bazer FW, Vallet JL, Roberts RM, Sharp DC, Thatcher WW. Role of conceptus secretory products in establishment of pregnancy, J Reprod Fertil. 76 (1986) 841-50. [4] Geisert RD, Yelich JV. Regulation of conceptus development and attachment in pigs, J Reprod Fertil. Page 15

ACCEPTED MANUSCRIPT 52 (1997) 133-49. [5] Vigano P, Mangioni S, Pompei F, Chiodo I. Maternal-conceptus cross talk--a review, Placenta. 24 (2003) S56-61. [6] Geisert RD, Brookbank JW, Roberts RM, Bazer FW. Establishment of pregnancy in the pig: II. Cellular

T

remodeling of the porcine blastocyst during elongation on day 12 of pregnancy, Biol Reprod. 27 (1982) 941-55.

IP

[7] Satterfield MC, Gao H, Li X, Wu G, Johnson GA, Spencer TE, et al. Select nutrients and their associated transporters are increased in the ovine uterus following early progesterone administration,

SC R

Biol Reprod. 82 (2010) 224-31.

[8] Bazer FW, Song G, Kim J, Dunlap KA, Satterfield MC, Johnson GA, et al. Uterine biology in pigs and sheep, J Anim Sci Biotechno. 3 (2012) 23.

[9] Bowen JA, Bazer FW, Burghardt RC. Spatial and temporal analyses of integrin and Muc-1 expression

NU

in porcine uterine epithelium and trophectoderm in vivo, Biol Reprod. 55 (1996) 1098-106. [10] Bowen JA, Bazer FW, Burghardt RC. Spatial and temporal analyses of integrin and Muc-1 expression in porcine uterine epithelium and trophectoderm in vitro, Biol Reprod. 56 (1997) 409-15.

MA

[11] Ka H, Jaeger LA, Johnson GA, Spencer TE, Bazer FW. Keratinocyte growth factor is up-regulated by estrogen in the porcine uterine endometrium and functions in trophectoderm cell proliferation and differentiation, Endocrinology. 142 (2001) 2303-10.

[12] Trowbridge JM, Rudisill JA, Ron D, Gallo RL. Dermatan sulfate binds and potentiates activity of

D

keratinocyte growth factor (FGF-7), J Biol Chem. 277 (2002) 42815-20.

TE

[13] Lennard SN, Gerstenberg C, Allen WR, Stewart F. Expression of epidermal growth factor and its receptor in equine placental tissues, J Reprod Fertil Suppl. 112 (1998) 49-57. [14] Wiley LM, Adamson ED, Tsark EC. Epidermal growth factor receptor function in early mammalian

CE P

development. BioEssays. 17(1995) 839-46. [15] Yoo HJ, Barlow DH, Mardon HJ. Temporal and spatial regulation of expression of heparin-binding epidermal growth factor-like growth factor in the human endometrium: a possible role in blastocyst implantation, Dev Genet. 21 (1997) 102-8.

AC

[16] Martin KL, Barlow DH, Sargent IL. Heparin-binding epidermal growth factor significantly improves human blastocyst development and hatching in serum-free medium. Hum Reprod. 13(1998) 1645-52. [17] Tamada H, Sakamoto M, Sakaguchi H, Inaba T, Sawada T. Evidence for the involvement of transforming growth factor-alpha in implantation in the rat, Life Sci. 60(1997) 1515-22. [18] Bai R, Kusama K, Sakurai T, Bai H, Wang C, Zhang J, et al. The Role of Endometrial Selectins and Their Ligands on Bovine Conceptus Attachment to the Uterine Epithelium During Peri-Implantation Period, Biol Reprod. 93 (2015) 46. [19] Lefevre F, Guillomot M, D'Andrea S, Battegay S, La Bonnardiere C. Interferon-delta: the first member of a novel type I interferon family, Biochimie. 80 (1998) 779-88. [20] Garlow JE, Ka H, Johnson GA, Burghardt RC, Jaeger LA, Bazer FW. Analysis of osteopontin at the maternal-placental interface in pigs, Biol Reprod. 66 (2002) 718-25. [21] White FJ, Ross JW, Joyce MM, Geisert RD, Burghardt RC, Johnson GA. Steroid regulation of cell specific secreted phosphoprotein 1 (osteopontin) expression in the pregnant porcine uterus, Biol Reprod. 73 (2005) 1294-301. [22] Vogiagis D, Salamonsen LA. Review: The role of leukaemia inhibitory factor in the establishment Page 16

ACCEPTED MANUSCRIPT of pregnancy, J Endocrinol. 160 (1999) 181-90. [23] Hirota Y, Daikoku T, Tranguch S, Xie H, Bradshaw HB, Dey SK. Uterine-specific p53 deficiency confers premature uterine senescence and promotes preterm birth in mice, J Clin Invest. 120 (2010) 803-15.

T

[24] Prins JR, Gomez-Lopez N, Robertson SA. Interleukin-6 in pregnancy and gestational disorders, J Reprod Immunol. 95 (2012) 1-14.

IP

[25] Dimitriadis E, Stoikos C, Tan YL, Salamonsen LA. Interleukin 11 signaling components signal transducer and activator of transcription 3 (STAT3) and suppressor of cytokine signaling 3 (SOCS3)

SC R

regulate human endometrial stromal cell differentiation, Endocrinology. 147 (2006) 3809-17. [26] Su L, Liu R, Cheng W, Zhu M, Li X, Zhao S, et al. Expression patterns of microRNAs in porcine endometrium and their potential roles in embryo implantation and placentation, PloS one. 9 (2014) e87867.

NU

[27] Chakrabarty A, Tranguch S, Daikoku T, Jensen K, Furneaux H, Dey SK. MicroRNA regulation of cyclooxygenase-2 during embryo implantation, Proc Natl Acad Sci U S A. 104(2007) 15144-9. [28] Kung JT, Colognori D, Lee JT. Long noncoding RNAs: past, present, and future, Genetics. 193(2013)

MA

651-69.

[29] Hu W, Alvarez-Dominguez JR, Lodish HF. Regulation of mammalian cell differentiation by long non-coding RNAs, EMBO Rep. 13 (2012) 971-83.

[30] Cesana M, Cacchiarelli D, Legnini I, Santini T, Sthandier O, Chinappi M, et al. A long noncoding

D

RNA controls muscle differentiation by functioning as a competing endogenous RNA, Cell. 147 (2011)

TE

358-69.

[31] Jalali S, Bhartiya D, Lalwani MK, Sivasubbu S, Scaria V. Systematic transcriptome wide analysis of lncRNA-miRNA interactions, PloS one. 8 (2013) e53823.

CE P

[32] Wang K, Long B, Zhou LY, Liu F, Zhou QY, Liu CY, et al. CARL lncRNA inhibits anoxia-induced mitochondrial fission and apoptosis in cardiomyocytes by impairing miR-539-dependent PHB2 downregulation, Nat Commun. 5 (2014) 3596. [33] Wang F, Ying HQ, He BS, Pan YQ, Deng QW, Sun HL, et al. Upregulated lncRNA-UCA1 contributes

AC

to progression of hepatocellular carcinoma through inhibition of miR-216b and activation of FGFR1/ERK signaling pathway, Oncotarget. 6 (2015) 7899-917. [34] Guttman M, Garber M, Levin JZ, Donaghey J, Robinson J, Adiconis X, et al. Ab initio reconstruction of cell type-specific transcriptomes in mouse reveals the conserved multi-exonic structure of lincRNAs, Nat Biotechnol. 28 (2010) 503-10. [35] Trapnell C, Williams BA, Pertea G, Mortazavi A, Kwan G, van Baren MJ, et al. Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation, Nat Biotechnol. 28 (2010) 511-5. [36] Sun L, Luo H, Bu D, Zhao G, Yu K, Zhang C, et al. Utilizing sequence intrinsic composition to classify protein-coding and long non-coding transcripts, Nucleic Acids Res. 41 (2013) e166. [37] Kong L, Zhang Y, Ye ZQ, Liu XQ, Zhao SQ, Wei L, et al. CPC: assess the protein-coding potential of transcripts using sequence features and support vector machine, Nucleic Acids Res. 35 (2007) W345-9. [38] Punta M, Coggill PC, Eberhardt RY, Mistry J, Tate J, Boursnell C, et al. The Pfam protein families database, Nucleic Acids Res. 40 (2012) D290-301. [39] Lin MF, Jungreis I, Kellis M. PhyloCSF: a comparative genomics method to distinguish protein Page 17

ACCEPTED MANUSCRIPT coding and non-coding regions, Bioinformatics. 27 (2011) i275-82. [40] Siepel A, Bejerano G, Pedersen JS, Hinrichs AS, Hou M, Rosenbloom K, et al. Evolutionarily conserved elements in vertebrate, insect, worm, and yeast genomes, Genome Res. 15 (2005) 1034-50. [41] Langmead B, Trapnell C, Pop M, Salzberg SL. Ultrafast and memory-efficient alignment of short

T

DNA sequences to the human genome, Genome biol. 10 (2009) R25. [42] Kozomara A, Griffiths-Jones S. miRBase: annotating high confidence microRNAs using deep

IP

sequencing data, Nucleic Acids Res. 42 (2014) D68-73.

[43] John B, Enright AJ, Aravin A, Tuschl T, Sander C, Marks DS. Human MicroRNA targets, PLoS Biol. 2

SC R

(2004) e363.

[44] Smoot ME, Ono K, Ruscheinski J, Wang PL, Ideker T. Cytoscape 2.8: new features for data integration and network visualization, Bioinformatics. 27 (2011) 431-2.

[45] Pennisi E. Long noncoding RNAs may alter chromosome's 3D structure, Science. 340(2013) 910.

NU

[46] Ponjavic J, Oliver PL, Lunter G, Ponting CP. Genomic and transcriptional co-localization of protein-coding and long non-coding RNA pairs in the developing brain, PLoS Genet. 5 (2009) e1000617.

MA

[47] Preker P, Nielsen J, Kammler S, Lykke-Andersen S, Christensen MS, Mapendano CK, et al. RNA exosome depletion reveals transcription upstream of active human promoters, Science. 322 (2008) 1851-4.

[48] Seila AC, Calabrese JM, Levine SS, Yeo GW, Rahl PB, Flynn RA, et al. Divergent transcription from

D

active promoters, Science. 322 (2008) 1849-51.

TE

[49] Sigova AA, Mullen AC, Molinie B, Gupta S, Orlando DA, Guenther MG, et al. Divergent transcription of long noncoding RNA/mRNA gene pairs in embryonic stem cells, Proc Natl Acad Sci U S A. 110 (2013) 2876-81.

CE P

[50] Core LJ, Waterfall JJ, Lis JT. Nascent RNA sequencing reveals widespread pausing and divergent initiation at human promoters, Science. 322 (2008) 1845-8. [51] Samborski A, Graf A, Krebs S, Kessler B, Reichenbach M, Reichenbach HD, et al. Transcriptome changes in the porcine endometrium during the preattachment phase, Biol Reprod. 89 (2013) 134.

AC

[52] Samborski A, Graf A, Krebs S, Kessler B, Bauersachs S. Deep sequencing of the porcine endometrial transcriptome on day 14 of pregnancy, Biol Reprod. 88 (2013) 84. [53] Lee SH, Zhao SH, Recknor JC, Nettleton D, Orley S, Kang SK, et al. Transcriptional profiling using a novel cDNA array identifies differential gene expression during porcine embryo elongation, Mol Reprod Dev. 71 (2005) 129-39. [54] Spencer TE, Bazer FW. Conceptus signals for establishment and maintenance of pregnancy, Reprod Biol Endocr. 2 (2004) 49. [55] Zhang H, Wang S, Liu M, Zhang A, Wu Z, Zhang Z, et al. Differential gene expression in the endometrium on gestation day 12 provides insight into sow prolificacy, BMC genomics. 14 (2013) 45. [56] Gu T, Zhu MJ, Schroyen M, Qu L, Nettleton D, Kuhar D, et al. Endometrial gene expression profiling in pregnant Meishan and Yorkshire pigs on day 12 of gestation, BMC genomics. 15 (2014) 156. [57] Kiewisz J, Krawczynski K, Lisowski P, Blitek A, Zwierzchowski L, Ziecik AJ, et al. Global gene expression profiling of porcine endometria on Days 12 and 16 of the estrous cycle and pregnancy, Theriogenology. 82 (2014) 897-909. Page 18

ACCEPTED MANUSCRIPT [58] Kim M, Seo H, Choi Y, Shim J, Kim H, Lee CK, et al. Microarray Analysis of Gene Expression in the Uterine Endometrium during the Implantation Period in Pigs, Asian-Australas J Anim Sci. 25 (2012) 1102-16. [59] Paiva P, Menkhorst E, Salamonsen L, Dimitriadis E. Leukemia inhibitory factor and interleukin-11:

T

critical regulators in the establishment of pregnancy, Cytokine Growth Factor Rev. 20 (2009) 319-28. [60] Giri JG, Newton RC, Horuk R. Identification of soluble interleukin-1 binding protein in cell-free

IP

supernatants. Evidence for soluble interleukin-1 receptor, J Biol Chem. 265 (1990) 17416-9. [61] Symons JA, Young PR, Duff GW. Soluble type II interleukin 1 (IL-1) receptor binds and blocks

SC R

processing of IL-1 beta precursor and loses affinity for IL-1 receptor antagonist, Proc Natl Acad Sci U S A. 92 (1995) 1714-8.

[62] Boucher A, Kharfi A, Al-Akoum M, Bossu P, Akoum A. Cycle-dependent expression of interleukin-1 receptor type II in the human endometrium, Biol Reprod. 65 (2001) 890-8.

NU

[63] Tsai SJ, Wu MH, Chen HM, Chuang PC, Wing LY. Fibroblast growth factor-9 is an endometrial stromal growth factor, Endocrinology. 143 (2002) 2715-21.

[64] Giudice LC, Milkowski DA, Fielder PJ, Irwin JC. Characterization of the steroid-dependence of

MA

insulin-like growth factor-binding protein-2 synthesis and mRNA expression in cultured human endometrial stromal cells, Hum Reprod. 6 (1991) 632-40. [65] Ostrup E, Bauersachs S, Blum H, Wolf E, Hyttel P. Differential endometrial gene expression in pregnant and nonpregnant sows, Biol Reprod. 83 (2010) 277-85.

D

[66] Klein C, Bauersachs S, Ulbrich SE, Einspanier R, Meyer HH, Schmidt SE, et al. Monozygotic twin

TE

model reveals novel embryo-induced transcriptome changes of bovine endometrium in the preattachment period, Biol Reprod. 74 (2006) 253-64. 373-422.

CE P

[67] Hoenderop JG, Nilius B, Bindels RJ. Calcium absorption across epithelia, Physiol Rev. 85 (2005) [68] Kligman D, Hilt DC. The S100 protein family, Trends Biochem Sci. 13 (1988) 437-43. [69] Swangchan-Uthai T, Lavender CR, Cheng Z, Fouladi-Nashta AA, Wathes DC. Time course of defense mechanisms in bovine endometrium in response to lipopolysaccharide, Biol Reprod. 87 (2012)

AC

135.

[70] Ka H, Seo H, Kim M, Choi Y, Lee CK. Identification of differentially expressed genes in the uterine endometrium on day 12 of the estrous cycle and pregnancy in pigs, Mol Reprod Dev. 76(2009) 75-84. [71] Geisert RD, Thatcher WW, Roberts RM, Bazer FW. Establishment of pregnancy in the pig: III. Endometrial secretory response to estradiol valerate administered on day 11 of the estrous cycle, Biol Reprod. 27 (1982) 957-65. [72] Ioannidis I, Dimo B, Karameris A, Vilaras G, Gakiopoulou H, Patsouris E, et al. Comparative study of the immunohistochemical expression of metalloproteinases 2, 7and 9between clearly invasive carcinomas and "in situ" trophoblast invasion, Neoplasma. 57(2010) 20-8. [73] Yanaihara A, Otsuka Y, Iwasaki S, Koide K, Aida T, Okai T. Comparison in gene expression of secretory human endometrium using laser microdissection, Reprod Biol Endocrinol. 2 (2004) 66. [74] Tapia A, Gangi LM, Zegers-Hochschild F, Balmaceda J, Pommer R, Trejo L, et al. Differences in the endometrial transcript profile during the receptive period between women who were refractory to implantation and those who achieved pregnancy, Hum Reprod. 23 (2008) 340-51. [75] Wu S, Sun H, Zhang Q, Jiang Y, Fang T, Cui I, et al. MicroRNA-132 promotes estradiol synthesis in Page 19

ACCEPTED MANUSCRIPT ovarian granulosa cells via translational repression of Nurr1, Reprod Biol Endocrinol. 13 (2015) 94. [76] Sang Q, Yao Z, Wang H, Feng R, Wang H, Zhao X, et al. Identification of microRNAs in human follicular fluid: characterization of microRNAs that govern steroidogenesis in vitro and are associated

SC R

IP

T

with polycystic ovary syndrome in vivo, J Clin Endocrinol Metab. 98 (2013) 3068-79.

Table 1 Number of reads and transcripts of RNA-seq Total reads

Clean bases

D12A D12B D12C D12K_A D12K_B D12K_C

85672786 87959926 84402788 90109252 83866148 91643288

10.7G 11G 10.56G 11.26G 10.48G 11.46G

Total mapped 84.28% 83.86% 83.04% 84.09% 82.55% 82.12%

AC

CE P

TE

D

MA

NU

Sample name

Page 20

Uniquely

mRNA

26279862 24813836 25497905 28045438 27353026 27512553

13448502(51.17%) 11322723(45.63%) 12565192(49.28%) 14416543(51.40%) 13757099(50.29%) 13323732(48.43%)

ACCEPTED MANUSCRIPT

AC

CE P

TE

D

MA

NU

SC R

IP

T

Table 2 Validation of RNA-seq results by using quantitative RT-PCR Mean RNA-seq Mean qPCR fold change Transcript ID correlation D12 D12K D12 D12K RNA-seq qPCR S100A8 10456.59 3.15 0.88 0.00 11.70 13.02 0.88 PRDX6 2120.46 683.51 0.46 0.12 1.63 1.98 0.80 STC1 3410.43 373.89 0.40 0.02 3.19 4.20 0.79 FGF9 76.63 9.87 0.52 0.08 2.96 2.78 0.80 PDK4 1893.38 782.73 0.87 0.16 1.27 2.40 0.78 CMTR1 1179.11 804.73 0.69 0.39 0.55 0.82 0.71 HEATR1 986.44 1415.26 1.79 4.14 -0.52 -1.21 0.76 TCONS_03397193 86.15 54.74 0.69 0.28 0.65 1.30 0.71 TCONS_00343128 70.81 4.77 6.17 3.57 3.89 0.79 0.75 TCONS_03397662 926.11 994.86 1.39 2.11 -0.10 -0.60 0.75 ssc-mir-34c 53460.00 67256.00 1.55 3.17 -0.33 -1.02 0.84

Page 21

ACCEPTED MANUSCRIPT

FEa

a

P-value

Genes

4.19 4.19

0.01 0.01

6 6

4.12

0.01

6

4.03 4.18 2.91 2.20 2.19

0.02 0.03 0.03 0.04 0.04

6 5 7 10 10

3.90

0.04

5

3.86 3.13

0.04 0.04

5 6

76.37

0.03

2

5.66 10.41 3.90 8.06 2.66

0.03 0.03 0.04 0.05 0.07

4 3 5 3 6

2.75

0.10

5

CR US

AC

CE P

TE D

MA N

Terms GOTERM_BP GO:0000280~nuclear division GO:0007067~mitosis GO:0000087~M phase of mitotic cell cycle GO:0048285~organelle fission GO:0008544~epidermis development GO:0000278~mitotic cell cycle GO:0007155~cell adhesion GO:0022610~biological adhesion GO:0031667~response to nutrient levels GO:0007398~ectoderm development GO:0051301~cell division GO:0008330~protein GOTERM_MF tyrosine/threonine phosphatase activity GO:0005179~hormone activity GO:0016836~hydro-lyase activity GO:0032403~protein complex binding KEGG_PATHWAY hsa02010:ABC transporters hsa04010:MAPK signaling pathway hsa04810:Regulation of actin cytoskeleton

IP

Category

T

Table 3 Significantly enriched GO terms and pathways from DAVID functional annotation

FE, fold enrichment Page 22

ACCEPTED MANUSCRIPT

T

Fig 1 The workflow to select lncRNA. We first assembled the transcripts using Cufflinks and scripture for each sample, then performed the

CR

IP

basic filtering by the five steps. At last, remove the potential coding transcripts by using the four software CPC, PfamScan, CNCI and PhyloCSF.

US

Fig. 2. Identification of non-coding lncRNAs by using CPC, PFAM, phyloCSF and CNCI. 1490 non-coding transcripts were selected using

MA N

four softwares that evaluate protein-coding transcripts and remove putative protein-coding transcripts. Fig. 3. Expression and relationship of lncRNAs with their neighboring coding genes. (a) Number of divergent and all-direction gene pairs formed by lncRNAs and their neighboring coding genes. Distribution of correlation of

TE D

each category of gene pairs (coding gene pair, lncRNA: coding gene pair and random pairs) in all neighbors (b) or in divergent neighbors (c), the left axis shows the density of gene pairs and the right axis shows the correlation of gene pairs. (d) Distribution of distance from one

CE P

transcription start sites (TSS) to another in all-direction lncRNA: coding gene pairs and coding gene pairs, and divergent direction lncRNA:

AC

coding gene pairs and coding gene pairs, the left axis shows the density of gene pairs and the right axis shows the distance to TSS. Fig. 4. GO analysis of differentially expressed lncRNAs and miRNAs. The differentially expressed genes (DEGs) are classified into three categories: cellular component (CC), molecular function (MF) and biological process (BP). The left axis shows the percentage of genes in a category, and the right axis shows the number of lncRNA (a) and miRNA (b) target genes.

Page 23

AC

CE P

TE D

MA N

US

CR

IP

T

ACCEPTED MANUSCRIPT

Fig. 1

Page 24

MA N

US

CR

IP

T

ACCEPTED MANUSCRIPT

AC

CE P

TE D

Fig. 2

Page 25

MA N

US

CR

IP

T

ACCEPTED MANUSCRIPT

AC

CE P

TE D

Fig. 3

Page 26

MA N

US

CR

IP

T

ACCEPTED MANUSCRIPT

AC

CE P

TE D

Fig. 4

Page 27

ACCEPTED MANUSCRIPT

T

Table 1 Number of reads and transcripts of RNA-seq

CR

IP

Table 2 Validation of RNA-seq results by using quantitative RT-PCR

US

Table 3 Significantly enriched GO terms and pathways from DAVID functional annotation

Highlights

Some genes in MAPK pathway may play a vital role in regulating embryo implantation.



XLOC_2604764 and XLOC_2604756 may be involved in regulating embryo implantation.



lncRNA-ssc-miR-132-mRNA may contributes to the establishment of maternal pregnancy.

AC

CE P

TE D

MA N



Page 28