Transcriptome profile of the human placenta

The human placenta is a particular organ that inseparably binds the mother and the fetus. The proper development and survival of the conceptus relies on the essential interplay between maternal and fetal factors involved in cooperation within the placenta. In our study, high-throughput sequencing (RNA-seq) was applied to analyze the global transcriptome of the human placenta during uncomplicated pregnancies. The RNA-seq was utilized to identify the global pattern of the gene expression in placentas (N = 4) from women in single and twin pregnancies. During analyses, we obtained 228,044 transcripts. More than 91% of them were multi-exon, and among them 134 were potentially unknown protein coding genes. Expression levels (FPKM) were estimated for 38,948 transcriptional active regions, and more than 3000 of genes were expressed with FPKM >20 in each sample. Additionally, all unannotated transcripts with estimated FPKM values were localized on the human genome. Highly covered splice junctions unannotated in the human genome (6497) were identified, and among them 30 were novel. To gain a better understanding of the biological implications, the assembled transcripts were annotated with gene ontology (GO) terms. Single nucleotide variants were predicted for the transcripts assigned to each analyzed GO category. Our results may be useful for establishing a general pattern of the gene expression in the human placenta. Characterizing placental transcriptome, which is crucial for a pregnancy’s outcome, can serve as a basis for identifying the mechanisms underlying physiological pregnancy, as well as may be useful for an early detection of the genomic defects.


Introduction
Next-generation sequencing (NGS) applied to the transcriptome (RNA-Seq) is a technology that uncovers broad amounts of novel information about the transcript (Lerch et al. 2012). RNA-Seq is a powerful tool suitable for a wide range of large-scale tissue-specific profiling, which includes analyses of transcriptome simultaneously for mutation, structural alteration, and expression, even the reconstruction of the entire transcriptome in the absence of a reference genome (Hui 2014). Despite the need for advanced bioinformatic and timeconsuming data analysis, the RNA-seq still provides a strong alternative to the use of microarrays in gene expression studies (Grada and Weinbrecht 2013). The discovery of specific genes crucial for placental development and function can help in understanding and identifying the mechanisms underlying both normal and pathological pregnancies. Detailing the placental tissue transcriptome could provide a valuable resource for genomic studies related to placental disease.
Electronic supplementary material The online version of this article (doi:10.1007/s10142-017-0555-y) contains supplementary material, which is available to authorized users.
Particularly in the context of established basis, the completion of the human genome sequencing project and the accessibility of the reference human genome enabled the application of genomic knowledge in clinical practice. The current progress of NGS technology is vertiginous, and genomic medicine driven by NGS may improve the diagnosis, prognosis, and therapy of many to-date incurable diseases (Ley et al. 2008;Dietz 2010).
Reproductive success is deeply linked with precise maternal-embryo communication, which results in correct placental functioning. The placenta is a primary organ essential for maintaining the proper development of the human fetus during pregnancy. Any disorder in this unique and delicate molecular crosstalk will have implications resulting in restricted progress in successful and proper embryo growth. Disturbances, if more severe, may lead to pregnancy failure, or even influence the developmental potential of offspring to adulthood (Fazeli and Pewsey 2008;Fan et al. 2016). Besides its hormonal and nutritive functions, the placenta also serves as a barrier against the maternal immune system, the gas exchange compartment, as well as the fetal waste disposal organ.
The nutritional status of the mother and infant from conception until early childhood, as well as maternal prenatal distress and compromise affect the early programming of the immune function, which might, in consequence, influence the neuro-development of the embryo (Marques et al. 2013).
Moreover, hormones produced by the placenta adjust maternal physiology to couple the body changes during the pregnancy's progress (Cross et al. 2003;Watson and Cross 2005). Additionally, significant differences in efficiency and growth, depending on the sex of the fetus, have been observed (Forsén et al. 1999;Eriksson et al. 2010). Differences in placental gene expression may also be correlated with the rate of growth and body length or weight as observed in males at birth, in comparison with females with equivalent placental size (Misra et al. 2009).
The properly functioning placenta depends on the finely tuned regulation of relevant genes that affect its structural integrity, and the proper growth and development of its structural components. Thus, any disturbances in gene expression may represent a major molecular mechanism underlying pathological pregnancies (Mikheev et al. 2008).
A profound penetration of the spatial and temporal changes of the placental transcriptome during development and pathology of the pregnancy should be considered essential in identifying the origins of diseases such as preeclampsia (PE) and intrauterine growth restriction (IUGR) that cause stillbirth or preterm birth. The consequences may extend into adult life, with a metabolic syndrome being a dramatic example. The advanced RNA-seq is an effective tool for analysis of the placental transcriptome to elucidate possible molecular mechanisms leading to these diseased states, and to identify candidate biomarkers (Cox et al. 2015).
The frequency of multiple gestations has increased over the last few decades, mainly due to the use of assisted reproductive technologies with the increased maternal age at childbearing. Twin pregnancies are at higher perinatal risk for both mothers and fetuses than singleton pregnancies (Chauhan et al. 2010;Hubinont et al. 2015).
The aim of the present study was to identify the transcriptome landscape of the human placenta during uncomplicated single and twin pregnancies to establish the pattern of normal placental gene expression for further comprehensive analyses.

Collection of placental samples and ethics statement
All clinical samples were collected at the Clinical Ward for Gynecology, Obstetrics and Oncological Gynecology at the Regional Specialist Hospital in Olsztyn following informed written consent from mothers. The experimental protocol was approved by the Bioethics Committee of the Warmia-Mazury Medical Chamber (OIL.164/15/Bioet) in Olsztyn, Poland.

Placental tissue collection
Placentas (N = 4) were collected from healthy women after uncomplicated pregnancy from term pregnancy (Table 1), who underwent a scheduled Cesarean section before the onset of labor. Samples of placental tissues were collected immediately after delivery; then, small pieces of placentas were frozen in liquid nitrogen. Preserved tissues were stored in the Laboratory of Molecular Diagnostics at −70°C until further analyses.
Total RNA extraction and purification RNA was isolated from the placental tissues using the Qiagen RNeasy Kit according to the manufacturer's recommendations. The RNase-Free DNase Set (Qiagen) was used to obtain high-quality RNA with DNA-ase digestion. RNA integrity (purity and concentration) was evaluated using the Tecan InfiniteM200 (Tecan Group AG). The RNA integrity number (RIN) was measured on a Bioanalyzer 2100 (Agilent). Only samples with the RIN >8 were subjected to further analysis. Fifteen microliters of total RNA extracted from each sample was delivered to OpenExome (Poland) in order to prepare a cDNA library and perform RNA-seq.

cDNA library construction and transcriptome sequencing
Double-stranded cDNA libraries were prepared from total extracted RNA using the TruSeq Stranded mRNA LT Sample Prep Kit v3 (Illumina) following manufacturer's instructions. The concentration of extracted RNA ranged from 10.07 to 24.06 ng/μl. Briefly, RNA samples were first purified with two oligo-dT selection (poly(A) enrichment using oligo-dT beds), and then fragmented and reverse transcribed into double-stranded complementary DNA. The resulting cDNA libraries of each sample were tagged with an indexed adapter prior to sequencing. The indexed libraries were diluted and pooled in equimolar ratios, and then the libraries were pair-end sequenced (101 cycles for read 1, 7 cycles for the index read, and 101 cycles for read 2) on the HiSeq2500 (Illumina), 2 × 100 bp reads were obtained.

Quality control and mapping to reference
The FastQC software (www.bioinformatics.babraham.ac. uk) was utilized for the evaluation of the quality control of the sequencing process. In order to remove Illumina adaptors and poly(A) stretches out of raw sequencing reads, the Trimmomatic tool (Bioinformatics, btu170) was used. Also, reads shorter than 50 bp or quality lower than Phred < 20 were removed from the dataset. Trimmed sequences were aligned to the human reference genome (Homo_sapiens.GRCh38.dna.primary_assembly. fa) with annotation (Homo_sapiens.GRCh38.80.gtf) file. Mapping was performed with the STAR (ver. 2.4) mapper with the following parameters: -outFilterType B y S J o u to u t F i l t e r M u l t i m a p N m a x 2 0 -alignSJoverhangMin 8 -alignSJDBoverhangMin 1o u t F i l t e r M i s m a t c h N m a x 9 9 9 -outFilterMismatchNoverLmax 0.04 -alignIntronMin 20 -alignIntronMax 1,000,000 -alignMatesGapMax 1,000,000. As a result, we obtained alignment of the trimmed reads to the reference genome (BAM files). The number of reads mapped to: exonic, intronic, UTRs, or i n t e r g e n i c r e g i o n s w a s q u a n t i f i e d u s i n g t h e CollectRnaSeqMetrics tool in the Picard 2.1.1 software (http://picard.sourceforge.net).

Transcriptional active regions
Clean reads were mapped on the human genome and as a result transcriptional active regions (TARs) were obtained. If TAR had an annotation (meaning assigned gene, non-coding region, processed transcript, miRNA, pseudogene, etc.), it was classified as annotated, if it did not have any annotation it was classified as unannotated.
Additionally, to specify the above analysis for each sample, cut-off FPKM (equal to 0.25) was applied to remove unannotated transcripts with low expression levels. Finally, 5320 TARs were expressed in more than one sample, 2742 TARs were expressed in more than two examined samples, and 880 unannotated TARs were expressed in all examined samples.

Alternative splicing events
Among all splice junctions (SJ.out.tab files) generated by STAR, novel exon-exon junction sites were selected. The number of uniquely mapped reads crossing the junction was set for 5×. RefGenome (bioconductor library in R) was used to overlap selected junctions with transcripts. Overlapped transcripts were merged with potential novel isoforms generated by Cufflinks-StringTie analyses. BLASTx was used to verify potential splice isoforms and indicate novel ones.

Gene ontology classification
The StringTie software was also used to generate additional annotations on the basis of mapped reads. Genes with FPKM >1 (sum of each sample) were selected, and the longest isoform of each transcript (unigene) was assigned to Gene Ontology (GO) terms in BLAST2GO software with BLASTx-fast; p value = 10e-5, human reference (Conesa et al. 2005). For those selected transcripts, InterproScan was applied to the predicted protein signature database including Pfam, Panther, and SMART (Hunter et al. 2009

Statistics of placenta transcriptome profile
The RNA-seq of the constructed cDNA libraries allowed the characterization of the transcriptome profile of the human placenta. High-throughput sequencing on an Illumina platform generated 2 × 119,560,140 total raw paired-end reads. After trimming, 2 × 109,363,183 reads with good quality were obtained. Among all clean reads, 86.76, 85.76, 84.29, and 87.72% were uniquely mapped for each sample, respectively. Moreover,~8.71% were mapped to multiple loci (2-20), 0.0125% to too many loci (>20), and only~0.0325% of reads were unmapped (Table 2).
For all input samples, the percentage distribution of aligned bases is mostly very similar:~3.4% of aligned bases were derived from intergenic regions,~34.85% of bases originated from UTR, and~55.8% from coding regions (Fig. 1). Significant discrepancies in aligned bases distribution were detected within intronic regions, between one and other three samples, 2 and~7.3% respectively (Fig. 1).

Transcriptional active regions
Cufflinks and StringTie generated 228,044 total transcripts, with more than 91% of them being multi-exon transcripts ( Table 3). The uniquely mapped reads were used to estimate the expression levels (FPKM) of 38,948 TARs (annotated and unannotated). To overview the gene expression, FPKM values were divided into appropriate intervals. Based on the FPKM value, more than 3000 of genes were expressed with FPKM >20 in each sample. The largest number of genes was expressed below 1 FPKM. The general pattern of the gene expression was similar in all samples (Fig. 2).
Among identified TARs, 29,514 were overlapping with the annotated transcripts from the NCBI RefSeq and Ensembl database. The remaining 9434 TARs were unannotated. The obtained RNA-seq data also enabled the prediction of 44,569 potentially novel isoforms, which on average is 5.86 isoforms per TAR (Table 3).
Based on a search against the human non-redundant (nr) database, unannotated regions were divided into four different groups of transcripts: over 90% similarity with nr database (558 transcripts), less than 90% similarity with ORF longer than 300 bp composed of less than 2 exons (754 transcripts), less than 90% similarity with ORF longer than 300 bp composed of at least 2 exons (138 transcripts), and less than 90% similarity with ORF shorter than 300 bp (8219 transcripts). Additionally, the last group, classified as non-coding RNA, was analyzed according to the length of sequences. We revealed 639.48 bp as the average length of transcripts longer than 200 bp, 2931 transcripts were long non-coding RNA (>500 bp), and the longest transcript was 11,781 bp in length (Fig. 3).
All of the unannotated transcripts with estimated FPKM values were localized on the human genome (Fig. 4).
Alternative splicing events STAR generated 239,458, 219,385, 235,672, and 281,038 splice junction sites separately for each sample, Hs_p14, Hs_p12, Hs_p9, and Hs_p3, respectively. Among all splice junctions, 6497 fulfilled the following conditions: novel exon-exon junction and uniquely mapped reads crossing the junction >5×. Overlapping with reference transcripts coordinates and merging with potential novel isoforms (StringTie) generated 2766 contigs. Among this group, BLASTx and non-redundant NCBI GenPept database confirmed 30 as novel splice junctions. Among them, pregnancy-associated plasma protein A (PAPPA) and hemoglobin subunit alpha (HBA) were chosen to broaden the analysis, namely, mapping and comparison with the human reference transcriptome and genome (Fig. 5).

Functional annotation based on GO categorization and KEGG annotation
To gain a better understanding of the biological implications, the assembled transcripts were annotated with gene ontology (GO) terms, based on the Bbest hit^BLASTx search in the Blast2GO pipeline. Unigenes with sum of FPKM >1 in all samples and FPKM >0.25 (in at least 3 samples) were annotated to three main GO categories: biological processes (82,751), cellular components (33,483), and molecular functions (32,625). Additionally, InterproScan predicted 3306 signatures in Panther, 16,994 in Pfam and 9145 in the SMART database (Online Resource 1). The majority of unigenes were assigned to biological processes GO category (51.3%), distributed in 27 subcategories, followed by cellular components (37%) and molecular functions (11.7%), with 19 and 17 subcategories, respectively.
The most overrepresented GO terms in biological processes were the cellular process (10,123) and single-organism process (8431). In the cellular component category, many transcripts were associated with intracellular membrane-bounded organelle (7424) and cytoplasm (7533). The molecular function category was mostly dominated by protein binding (7419), ion binding (4305), and organic cyclic compound binding (4020).
The top 10 represented GO terms ( Fig. 6a and Online Resource 2) for the biological process were transcription (1211), regulation of transcription (799), and positive regulation of transcription from RNA polymerase II promoter (618). For molecular function, the top represented terms ( Fig. 6b and Online Resource 2) were protein binding (2427), metal ion binding (1329), and ATP binding (1139). Lastly, the top cellular component GO terms ( Fig. 6c and Online Resource 2) were cytosol (1960), extracellular exosome (1746), and nucleoplasm (1388). Due to our concern about widely understood reproduction, we choose to present GO categories associated with this process. Thorough GO analyses revealed that 53 genes were a s s i g n e d t o t h e em b r y o d ev el op men t c ate go ry. Furthermore, 29 genes were annotated to embryo implantation, 26 genes were categorized to placenta development, and 17 genes from this category were annotated to embryonic placenta development, whereas 9 genes were annotated to maternal placenta development (Fig. 7).
Additionally, molecular pathways were annotated from the Kyoto Encyclopedia of Genes and Genomes database (KEGG). Purine metabolism, biosynthesis of antibiotics, and pyrimidine metabolism were the most numerous, with 260, 219, and 120 sequences in the pathway, respectively (Online Resource 3). Variant (insertion, deletion, and substitution) detection was conducted for the contigs exclusively found for each GO term: embryo development, system development, embryo implantation, and embryonic and maternal placenta development (Table 4).

Discussion
In the current study, we performed RNA-Seq to characterize the transcriptome profile of the human placenta term pregnancies ended by cesarean section. Our thorough analyses of placental transcriptome permitted the identification of transcripts' profile, novel alternative junction events, and potential SNV within the entire placental transcriptome. Results obtained in   TAR transcriptional active regions our study complement current findings in the field of molecular profiling of pregnancy-depended transcriptome changes in various tissues (Adams et al. 2011;Sharp et al. 2016). We obtained 119,560,140 total raw reads with 196-97 bp average read length from four placental term tissue, whereas Kim et al. (2012) generated paired-end data with 54 and 72 bp in read length and 23-33 million reads of each lane, for a total of 50-60 million paired-end reads per tissue for each of the placental tissues (amnion, chorion, and decidua). We obtained a high mapping rate with uniquely mapped reads in the range of 84.29-87.72% of reads aligned to the reference genome and only 1.53, 6.95, 7.44%, and 4.67% (for each of sample) of total unmapped reads, which indicates successful cDNA library construction and good sequencing quality of the human placental transcriptome. Using the uniquely mapped read pairs,  estimated the expression levels of 22,523 protein-coding genes for single placental tissue, and approximately half of the genes are expressed with FPKM >1. Among total expressed TARs, our research revealed 14,785-16,202 protein-coding genes with levels of expression FPKM >1 for four placental tissues. Our data show that the largest number of genes was expressed below 1 FPKM (Fig. 3).
In our research, we calculated FPKM values for every transcript from each sample to generate an overview on placental transcriptome. The highest FPKM values were determined for high temperature requirement protease A1 (HTRA1), endothelial PAS domain protein 1 (EPAS1), and Pleckstrin homologylike domain, family a, member 2 (PHLDA2).
The physiological function of HTRA1 remains unclear. However, human HTRA1 has been linked with the pathogenesis of several abnormalities, such as osteoarthritis, rheumatoid arthritis, and even cancers (Hasan et al. 2015). Additionally, HTRA1 concentration increases in the maternal plasma of gestations complicated by preeclampsia with intrauterine growth restriction compared with control (D'Elia et al. 2012). EPAS1 encodes a transcription factor involved in the response to fluctuating oxygen levels. There are conflicting reports concerning associations of EPAS1 mutations and expression level with abnormalities during pregnancy. However, EPAS1 is considered a candidate gene as a potential biomarker of the susceptibility, early detection, and/or individualized maternal-infant care in case of preeclampsia (Founds et al. 2011). PHLDA2 regulates fetal growth, placental development, and placental lactogen production (Tunster et al. 2013;Jensen et al. 2014;Janssen et al. 2016). Abnormally increased placental PHLDA2 expression is associated with fetal growth restriction in human pregnancies and/or low birth weight (Jensen et al. 2014). Increased PHLDA2 expression is also reported in cases of spontaneous miscarriage or fetal death (Doria et al. 2010). Moreover, increased placental PHLDA2 expression correlates with reduced fetal movements (RFM) pregnancies. Because of its impact on pregnancy outcome, PHLDA2 is postulated as a placental prenatal diagnostic biomarker identifying RFM pregnancies (Janssen et al. 2016). For our study, we selected patients with uncomplicated pregnancies, and thus we cannot exclude that high level of placental mRNA for HTRA1 does not affect its serum level. However, we also cannot exclude that some undiagnosed or unnoticeable abnormality might have occurred in samples with high expression of HTRA1, PHLDA2, and EPAS1.
In our study, searching unannotated regions against the nr database revealed 134 transcripts with ORF longer than 300 bp composed of at least 2 exons and similarity less than 90%. In similar studies, 604, 1007, and 896 novel TARs in amnion, chorion, and decidua, respectively, not overlapping with the annotated transcripts from the NCBI RefSeq, UCSC, Ensembl, and Vega database were identified ). Thus, taking into consideration the spatio-temporal changes of the transcriptome landscape, the differences between the obtained results and previous studies are inescapable. Fig. 3 Distribution of noncoding RNA according to length (bp) in the human placenta Fig. 4 Comparison of unannotated transcripts' distribution and expression level across the human genome. Four outer tracks depict the gene expression level (FPKM ranges from 0 to 20) of each placental sample. The most outer -Hs_p14 followed by Hs_p12, Hs_p9, and Hs_ p3. Colored tags represent four different groups of transcripts: >90% similarity with the Human non-redundant database (gray circles), <90% similarity with ORF <300 bp (green circles), <90% similarity with ORF >300 bp composed of <2 exons (purple circles), <90% similarity with ORF >300 bp composed of >2 exons (red triangles). The innermost circle shows the percentage of every group to the total number of unannotated transcripts (9669) Additionally, our studies uncovered 6497 splice junctions with high-quality coverage, and among this group 30 were verified as novel splice junctions.
Among the identified 30 novel splice junctions, we considered pregnancy-associated plasma protein A (PAPPA) and hemoglobin subunit alpha (HBA) the most interesting in the context of our research. PAPPA is a metalloprotease belonging to the metzincin superfamily of zinc peptidase (Leguy et al. 2014). It is an active homodimer (dPAPPA), circulating at very low levels in non-pregnant women and men, and regulates cell differentiation and proliferation by cleaving insulin-like growth factor binding proteins 4 and 5 (D'Elia et al. 2012;Coskun et al. 2013). During pregnancy, PAPPA is produced at high levels by the placenta and circulates as a heterotetrameric complex (htPAPPA). Decreased levels of this complex are associated with adverse pregnancy outcomes such as intrauterine growth restriction (IUGR), preterm delivery, miscarriage, preeclampsia, or fetal aneuploidy (Kirkegaard et al. 2010;Ranta et al. 2011). According to The World Health Organization, deletions of the α-globin genes are the most common genetic disorders in the world, and afflict 20% of the world's population (Modell and Darlison 2008). The deletion of one or two α-globin genes may provide positive effects such as protection against malaria (Harteveld and Higgs 2010). However, the deletion of three genes causes thalassemia, which often results in significant anemia. Moreover, the deletion of all four α-globin genes results in fetal death and significant risk of maternal morbidity or mortality (Leung et al. 2008). The development of high-throughput RNA-seq and bioinformatic analyses enabled the detection of alternative splicing events (Wang et al. 2009). In general, genes contain multiple introns, with an average of 7.8-9.0 introns per gene in vertebrates. Previous RNA-seq analyses revealed that >80% of mammalian genes include introns, and >95% of transcripts in humans are alternatively spliced (Pan et al. 2008;Wang et al. 2008). Alternative splicing can generate various transcripts and thus multiple protein isoforms from a single gene. This mechanism allows the modulation of gene function, phenotype, and eventually cause diseases (Mourier and Jeffares 2003). Transcriptome and proteome diversity can be generated by exon skipping, alternative 5′ donor sites, and intron retention (Li et al. 2014). Regulation of alternative splicing (AS) events in multi-exon genes not only affects expression level or protein function, but can also lead to differences in AS between human individuals (Wang et al. 2008;Lu et al. 2012). Changes in coding sequence may lead to disturbances in functions, due to the presence or absence of specific protein domains. Alternative splicing events may also be present in promoter and UTR regions affecting expression level and isoform localization, respectively . Alternative splicing events regulate gene expression, which is essential for cell differentiation and organism development (Horiuchi and Aigaki 2006;Pan et al. 2008;Nilsen and Graveley 2010). It is mostly unknown how splice variants differ in terms of properties, but it is still incontestable that alternative splicing may increase the functional diversity in positive and negative manners (Polo-Parada et al. 2004;An et al. 2008). Alternative splicing as a universal regulatory mechanism has a significant impact on the tissue-specific transcriptome landscape (Wang et al. 2008). The example of the placental matador gene splice isoform that is responsible for trophoblast cell death in preeclampsia (Soleymanlou et al. 2005) shows how placental transcriptome may affect embryo development. Additionally, this underlines the importance of research aimed at characterizing placental transcriptome. We can assume that the novel placental transcripts identified in this study (such as PAPPA and HBA) may contribute to pregnancy course, but further functional analyses are required to state how they affect pregnancy outcome. To gain more insight into placental functioning during uncomplicated pregnancy, we performed functional annotation analysis of placental transcripts that may be useful in understanding the functional significance of the placental transcripts, and in the future may constitute the answer to how they are implicated in placental biology and/or pregnancy disorders. In a similar study, differentially expressed placental enriched transcripts were mostly assigned to the following GO categories: multicellular organismal process, singlemulticellular organism process, and response to stress (Hughes et al. 2015). We performed GO analyses for all placental transcripts with FPKM >1 to generate an overview of the placental landscape. In our analyses, transcripts annotated to biological processes were distributed in 16 subcategories. Over the represented category was the cellular process. However, most enriched GO terms (Hughes et al. 2015) were also represented in our results. In the general overview, it may be observed that placental expression is very diversified and numerous, and distinct GO categories were represented (Figs. 6 and 7). This is probably the effect of wide placental  functions, including the transport of gases, nutrients and waste products, hormone production, and protection of the fetus from the maternal immune system (Kay et al. 2011). Next-generation sequencing, including deep sequencing of the complementary RNA-seq, allows for global highthroughput analyses, leading to the identification of SNVs within known genes, novel transcripts, alternative splicing variants, trans-splicing events, and to adjusting predicted gene models (Wang et al. 2008(Wang et al. , 2009. The obtained results regarding SNVs within genes associated with reproduction may constitute a basis for further genome-wide association (GWA) studies (Appels et al. 2013;Akpinar et al. 2017). The relationship between a specific genotype and a phenotype can be used to predict genes that may be correlated with observable traits in animals. GWA studies based on NGS results do not depend upon prior knowledge of the biological pathways, and this provides the opportunity to discover SNPs associated with phenotype (Liaoa and Leeb 2010). However, in the case of our results, it is possible to predict SNVs that may be involved in proper placental development or that induce pathology during pregnancy, which can seriously affect pregnancy outcome. In the light of recent research showing that~20% of infertility has an unknown background (Iammarrone et al. 2003), our study is in line with the search for genetic markers that may be associated with idiopathic infertility.
In conclusion, we characterized the transcriptomic profile of placental tissue and identified novel transcripts expressed in term human placenta that may be involved in pregnancy course. These results may provide a basis for further analyses of the pathways essential for proper placental development and functioning that contribute to the pregnancy outcome. Overall, these results provide new resources for the identification of mechanisms regulating placental physiology during non-complicated pregnancy.
Compliance with ethical standards Data availability statement The sequencing data from this study have been submitted (http://www.ncbi.nlm.nih.gov/sra) to the NCBI Sequence Read Archive (SRA) under accession no. BioProject ID: PRJNA326064. As well, all identified cDNA sequences of the novel transcripts have been deposited in the GenBank (KX533474-KX533503).
Funding This study was supported by the Department of Human Physiology, Faculty of Medical Science (25.610.001-300), University of Warmia and Mazury in Olsztyn.

Conflict of interest
The authors declare that they have no conflict of interest.
Ethical approval All procedures performed in studies involving human tissues were in accordance with the ethical standards of the institutional committee and with the 1964 Helsinki declaration and its later amendments or comparable ethical standards.
Open Access This article is distributed under the terms of the Creative Comm ons Attribution 4.0 International License (http:// creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.