RNA sequencing of chorionic villi from recurrent pregnancy loss patients reveals impaired function of basic nuclear and cellular machinery

Recurrent pregnancy loss (RPL) concerns ~3% of couples aiming at childbirth. In the current study, transcriptomes and miRNomes of 1st trimester placental chorionic villi were analysed for 2 RPL cases (≥6 miscarriages) and normal, but electively terminated pregnancies (ETP; n = 8). Sequencing was performed on Illumina HiSeq 2000 platform. Differential expression analyses detected 51 (27%) transcripts with increased and 138 (73%) with decreased expression in RPL compared to ETP (DESeq: FDR P < 0.1 and DESeq2: <0.05). RPL samples had substantially decreased transcript levels of histones, regulatory RNAs and genes involved in telomere, spliceosome, ribosomal, mitochondrial and intra-cellular signalling functions. Downregulated expression of HIST1H1B and HIST1H4A (Wilcoxon test, fc≤0.372, P≤9.37 × 10−4) was validated in an extended sample by quantitative PCR (RPL, n = 14; ETP, n = 24). Several upregulated genes are linked to placental function and pregnancy complications: ATF4, C3, PHLDA2, GPX4, ICAM1, SLC16A2. Analysis of the miRNA-Seq dataset identified no large disturbances in RPL samples. Notably, nearly 2/3 of differentially expressed genes have binding sites for E2F transcription factors, coordinating mammalian endocycle and placental development. For a conceptus destined to miscarriage, the E2F TF-family represents a potential key coordinator in reprogramming the placental genome towards gradually stopping the maintenance of basic nuclear and cellular functions.


Elective (surgical) termination of pregnancy (ETP group):
The subjects (n=8 in the discovery RNA-Seq study; n=24 for TaqMan RT-qPCR experiment) had experienced normal first trimester pregnancies with no maternal or fetal clinical complications until the elective termination of the pregnancy (mean maternal age 25.5±5.9 years, and 27.3±6.4 respectively; gestational age 64.6±13.9 days and 64.6±11.6 days respectively; Table 1 and Supplementary Table S1). None of the patients had experienced any clinically confirmed pregnancy losses in their reproductive history. The RPL cases recruited to the placental RNA-Seq/miRNA-seq discovery study had experienced five (RPL1) and six (RPL2) clinically confirmed pregnancy losses before the index case. The RPL placental samples were collected at surgical removal of conceptus shortly (within 24 h) after detection of fetal death. In both cases, the most recent ultrasound scan confirming the heart beats and normal growth of the fetus had been performed five (RPL1) and seven (RPL2) days before the event. In both patients diagnostic hysteroscopy, performed after index pregnancy revealed partial (up to 1 cm) uterine septum that was removed during the procedure. The anomaly has not been detected at multiple U/S performed during previous years before the index pregnancy. After the index pregnancy and removal of uterine septum, both patients have had one term delivery and further pregnancy losses. Taqman RT-qPCR experiments included additional placental samples (ETP n=24; RPL n=12). For these samples, the washing step in PBS was followed by immediate visual separation of chorionic villi, which were snap-frozen in liquid nitrogen or placed into RNAlater solution (Ambion Inc, Life Technologies) and kept at -80°C until RNA isolation [1].

Placental samples of normal term pregnancies
The RNA-Seq dataset representing term placental transcriptome in normal, uncomplicated pregnancies (Term norm) was derived from previous published report [2] and utilized as a reference to compare with the 1 st trimester placental gene expression. miRNA-seq data of the identical samples (n=8) was generated in the current study.
Term placental samples representing normal pregnancy originate from the REPROgrammed fetal and/or maternal METAbolism (REPROMETA) sample collection [2][3][4]. The study group in the current study (n=8, Table 1) was initially reported in our recent study comprised of uncomplicated term pregnancies (maternal age 29.3±7.9, gestational age 278.6±11.5), which resulted in the delivery of a newborn with normal birth weight for its gestational age (birth-weight between 10 th -90 th centile, [5]).
Cases with documented fetal anomalies, chromosomal abnormalities, families with history of inherited diseases and patients with known pre-existing diabetes mellitus, chronic hypertension and chronic renal disease were excluded. Placentas (stored at +4°C) were sampled within 1 h after vaginal delivery (n=5) or elective caesarean section without labor (n=3). Full-thickness blocks of 2-3 cm were taken from the middle region of the placenta. Collected tissue samples were washed with 1x PBS to remove contamination of maternal blood, placed immediately into RNAlater solution (Ambion Inc, Life Technologies) and kept at -80°C until RNA isolation.

RNA extraction and purification
Total RNA was extracted from 200-300 mg of homogenized placental tissue using TRIzol reagent (Invitrogen, Life Technologies). For RNA-Seq protocol, it was further purified with RNeasy MinElute columns (Qiagen, Netherlands) and for Taqman RT-qPCR applications with NucleoSpin ® II Isolation Kit (Macherey-Nagel GmbH & Co. KG, Düren, Germany) according to the manufacturers' protocols. Purity level and concentration of isolated total RNA was measured using NanoDrop ® ND-1000 UV-Vis spectrophotometer (Thermo Fisher Scientific Inc., USA) and RIN (RNA integrity number) was estimated using Agilent 2100 Bioanalyzer (Agilent Technologies, USA). control; reads were filtered for adaptor, rRNA and mtDNA sequences as well as homopolymer stretches using custom python scripts; read alignment was performed with TopHat version 2.0.3 [6] using bowtie version 0.12.7 [7]; transcript quantification was conducted with Cufflinks v 2.0.2 [8] with reference annotation (measured as FPKM) and gene expression was quantified by htseq-count [9] (as raw read counts). Human genome assembly (GRCh37.p7/hg19) from Ensembl v67 was used as a reference. Read alignment to genomic regions was estimated with Picard v1.63 (http://picard.sourceforge.net).

Library preparation and basic bioinformatics of raw data for miRNA-seq:
Initial small-RNA libraries were prepared from 1 µg total RNA (TruSeq Small RNA kit, Illumina), followed by microRNA enrichment (Caliper LabChipXT, PerkinElmer) from the 147 bp+/-5% PCR amplification fraction according to manufacturer's protocols. Libraries were sequenced on Illumina HiSeq2000.
Initial analysis of microRNA sequencing data included trimming of reads by removing the 3' adaptor sequence and discarding reads with the trimmed length <14 bp, mapping the remaining reads with bowtie (v0.12.7) to the reference genome (GRCh37.p7/hg19) allowing a maximum of 1 mismatch. Calculation of the raw read counts and normalized RPKM was performed using a modified version of Emir bioinformatic pipeline [10], allowing reads that map to genome multiple times (maximum 20) and using mature miRNA expression database (mirBASE v18).

Statistical analysis
Differential expression in RNA-Seq and miRNA-seq data was tested using DESeq [11] and DESeq2 [12] packages for R [13]. Read counts from htseq-count were used as input (number of Ensembl v67 genes with read counts n = 53,893). Genes with mean normalized expression < 50 reads in all samples (n = 39,109 RNA-Seq; n = 555 miRNA-seq) were considered as transcriptional noise and filtered out from the analysis. After the exclusion of genes with negligible placental transcript counts, 14,784 genes from RNA-Seq and 312 miRNAs entered differential expression testing. Built-in normalization algorithms of DESeq and DESeq2 were used. Outlier detection and handling was performed using the default method in DESeq. In DESeq2 outliers were replaced using the replaceOutliersWithTrimmedMean function with default Cook's distance cutoff. Statistical testing indicated that the two software packages, DESeq and DESeq2 differ substantially for their sensitivity in assessment of differential expression. Compared to the seminal DESeq package, analysis with the more recently developed DESeq2 programme produced a markedly higher number of significant results for all conducted differential expression tests with our data.
Thus, in the current study a more stringent level of significance was imposed on the test results of DESeq2. A gene was considered as differentially expressed, when the statistical tests simultaneously satisfied the following thresholds: FDR<0.1 for DESeq and FDR<0.05 for DESeq2.

Taqman RT-qPCR experiments and statistical analysis of RT-qPCR data
For an independent experimental validation, gene expression of five protein-encoding genes was analyzed using Taqman RT-qPCR gene expression assays (Applied Biosystems, Life Technologies; Supplementary  Gene expression was quantitated by biplex RT-qPCR of the target gene and housekeeping gene YHWAZ sequence using pre-made TaqMan Gene Expression Assays (Hs03044281_g1, Applied Biosystems, Life Technologies). YHWAZ was selected as the most stable ready-to-use housekeeping gene based on data from RNA-Seq and literature sources [14,15]. cDNA was synthesized from 1 µg total RNA according to the manufacturer's instructions (SuperScript III Reverse Transcriptase, Life Technologies).
All qPCR reactions were performed in triplicate in 384 micro-well plates in ABI 7900HT Real-