Comparison of the myometrial transcriptome from singleton and twin pregnancies by RNA-Seq

Preterm birth is recognized as the primary cause of infant mortality worldwide. Twin pregnancies are significantly more at risk of preterm birth than singleton pregnancies. A greater understanding of why this is and better modes of treatment and prevention are needed. Key to this is determining the differing pathophysiological mechanisms of preterm birth in twins, including the role of the myometrium and premature uterine contraction. We performed RNA sequencing (RNA-Seq) of human myometrium from singleton and twin pregnancies at term (> 37+0 weeks) and preterm (< 37+0 weeks), collected during pre-labour Caesarean Section. RNA-Seq libraries were prepared from polyA-selected RNA and sequenced on the Illumina HiSeq 4000 platform. Differentially expressed genes (DEGs), GO (Gene Ontology) and KEGG (Kyoto Encyclopedia of Genes and Genomes) pathway enrichment were conducted using R software. Significance was determined with a false discovery rate–adjusted P value of <0.05. Only 3 DEGs were identified between gestational age-matched singleton and twin myometrium and only 1 DEG identified between singleton term and twin preterm tissues. Comparison of singleton preterm myometrium with twin term myometrium however, revealed 75 down-regulated and 24 up-regulated genes in twin myometrium. This included genes associated with inflammation and immune response, T cell maturation and differentiation and steroid biosynthesis. GO and KEGG enrichment analyses for biologically relevant processes and functions also revealed several terms related to inflammation and immune response, as well as cytokine-cytokine receptor interaction and chemokine receptor signalling. Data indicate that little or no differences exist in the transcriptome of singleton and twin myometrium when matched for gestational age. The significant up- and down-regulation of genes identified between preterm singleton and twin myometrium at term may point to transcriptome changes associated with the chronic levels of uterine stretch in twin pregnancy or genes associated with the myometrium transitioning to labour onset.

Introduction Preterm birth is defined as delivery before 37 weeks gestation and occurs at a rate of around 11% of all births worldwide [1]. It is the single biggest cause of neonatal mortality and morbidity with almost half of preterm births being preceded by preterm labour (PTL). Twin pregnancies represent approximately 3% of all births yet, in contrast to singleton pregnancy, twin pregnancy is characterised by much higher rates of preterm birth: Nearly 60% of all twins are born before 37 weeks and around 20% deliver before 34 weeks [1,2]. Twin pregnancies therefore contribute disproportionately to the rates of preterm birth and adverse neonatal and infant outcomes and place a greater public health, economic and emotional burden on society.
As with singleton pregnancies, the aetiology of preterm birth in twins is likely to be multifactorial but risk factors predisposing to preterm birth in twins may be different to that of a singleton pregnancy [3,4]. Infection or short cervix are known risk factors for preterm birth [5] and are likely to be common to both pregnancies, whilst others such as uterine (over) distension and early activation of contractile pathways may be more significant in twins. Endocrine differences may also exist, resulting in part from double the quantity of fetal hormones from having more than one fetus, each producing hormones in the fetal hypothalamic-pituitary-adrenal-placental axis and in the case of dichorionic twins, the effects of increased overall placental mass [3].
Uterine over-distention or 'stretch' which refers to the mechanical stress on the uterine wall following the increase in uterine volume, has long been thought to contribute to preterm labour onset: In humans, placement and inflation of a balloon extra-amniotically can induce abortion in the second trimester and induce uterine contractions at term or postpartum [6][7][8]. Interestingly, women randomized to larger volumes of balloon have significantly shorter time from induction to delivery [9]. Increased intrauterine volume using intraamniotic balloons in non-human primates also initiated preterm labour in 50% of animals [10]. Hence in twin pregnancy, or pregnancies complicated by polyhydramnios (excess amniotic fluid) or macrosomic (large) fetus, the resulting over-distention of the myometrium may contribute to the higher rates of PTL in these patients. Stretching of human myometrium in vitro is also linked to upregulation of contraction associated proteins such as cyclooxygenase 2 (COX-2) [11] and the oxytocin receptor (OTR) [12], pro-inflammatory factors such as IL-8 [13] as well as release of cytokines from myometrium [14] which are key inflammatory mediators involved in labour onset (preterm or term).
As is the case for singleton pregnancy, predicting and preventing spontaneous PTL and birth are important goals for optimising outcomes in twin pregnancies. Once labour has been initiated, there are a limited number of effective and safe treatments (tocolytics) to prevent birth [15]. Current treatments however, do not substantially improve neonatal and maternal outcomes [16]. Preterm labour prevention strategies such as progesterone for women with short cervix and history of PTL have shown some promise in singleton pregnancy [17]. However, all interventions appear to be less effective in twin pregnancy [18]. There is therefore a need for improvement in current tocolytic treatments or a targeted/personalised treatment approach for singleton and twin pregnancy groups. However, this requires a greater understanding of the different mechanisms (if any) which bring about PTL in twins and singletons.
In myometrial function studies, small but significant differences in contractility of twin myometrium have been observed including a heightened response to oxytocin [19] and differences in the response of twin pregnancy myometrium to some tocolytics and progesterone when contractions were stimulated with oxytocin [20,21]. Together this suggests that there may be differences in the expression pattern of a number of proteins associated with contraction and its modulation in twins. Recently, TREK1 expression was implicated in maintaining uterine quiescence following prolonged stretch in twin pregnancy myometrium [22]. Studies involving gene expression in twin pregnancy myometrium however are few in number [10,22,23] and the mechanisms of PTL in twins is relatively understudied.
Previous studies examining the expression of mRNA or proteins in both humans and animals, show that there are marked changes in myometrial gene expression which precede the onset of labour [24]. A number of studies have sought to explore these changes at the whole genome level using DNA microarrays in both humans [24][25][26][27][28][29][30] and rodents [31][32][33], or more recently have adopted an RNA sequencing (RNA-Seq) approach [34][35][36]. RNA-Seq provides an open, unbiased method to explore transcriptomic changes at the genome level for a given cell or tissue such as in response to treatment or differences between patient groups. To date in human myometrium, it has been successfully applied to explore changes in the transcriptome in singleton pregnancies with labour onset [34,36]. No study has yet compared the gene expression profiles of myometrium between singleton and twin pregnancies.
In this paper we present the results of the first RNA-Seq study comparing gene expression in myometrium from pre-labour twin and singleton pregnancies matched for preterm or term gestation. Our objective was to identify differences in the myometrial transcriptome between singleton and twin pregnancies which would point to a unique gene activation pattern (possibly attributed to increased uterine stretch) and provide insight into the early activation of contraction pathways that may underpin aberrant contractility in preterm labour in multiple pregnancy.

Myometrial sample collection
All samples were obtained with informed written consent and the study was approved by the local research ethics committee (North West 3 Research Ethics Committee-Liverpool East, Ref. 10/H1002/49) and had Institutional review board approval from the Research and Development department at Liverpool Women's Hospital NHS Foundation Trust. Women were approached by a research midwife during Caesarean section (CS) booking appointments in pre-operative clinic or by a doctor during a routine scan appointment in the fetal medicine unit at Liverpool Women's Hospital. Recruitment took place between February 2014 and October 2017.
Criteria for inclusion in the study was singleton or twin pregnancy booked for elective (prelabour) CS delivery between 32 and 39 weeks of gestation where the estimated due date had been determined by ultrasound scan during pregnancy booking (10-12weeks). Women were excluded if they were found to have any medical condition e.g. diabetes, hypertension or preeclampsia or receiving any medication.
Samples were obtained by knife biopsy from the upper lip of the lower uterine segment during surgery, following delivery of the baby/ies and prior to oxytocin administration. Samples were immediately placed into Hanks balanced salt solution, transferred to the laboratory where they were inspected under a dissection microscope. Biopsies were cleared of small blood vessels, decidua, membrane and scar tissue (if present) to ensure a homogeneous cell type i.e. myometrium, was collected. Small pieces of myometrium were then placed into RNA later for 24 hours at 4˚C and stored at -80˚C prior to RNA extraction.
Samples were divided into 4 groups (n = 6 per group) according to gestation at delivery and multiplicity of pregnancy: Singleton Preterm (SP) and Twin Preterm (TP) delivering before 37 +0 weeks gestation and Singleton Term (ST) and Twin Term (TT) delivering after 37+0 weeks gestation. Reasons for CS included: Previous CS or previous traumatic vaginal delivery (PTVD, n = 5), breech presentation (n = 4), maternal request (n = 7), poor fetal growth (n = 3), previous preterm delivery resulting in CS delivery (n = 2) or maternal reason unrelated to pregnancy (n = 3). A summary of maternal characteristics according to gestational age group are presented in Table 1.
For both preterm and term groups, there was no difference in mean/median values of maternal age, maternal BMI, parity or gestational age at delivery between singleton and twin pregnancy women, indicating that the groups were well-matched demographically. As expected however, neonatal birthweight (which was combined birthweight in the case of twins), was significantly greater in both twin pregnancy groups compared to their gestational age-matched singleton counterparts. Overall in comparing all four groups, neonatal birthweight in Twin Term > Twin Preterm > Singleton Term > Singleton Preterm. The indications for CS delivery were also significantly different between groups. Preterm groups tended to require CS delivery for issues related to the fetus or maternal medical conditions (not related to pregnancy), whilst term deliveries were for previous CS, PTVD, breech presentation or maternal request (Table 1).

RNA extraction and quality assessment
Eighty to one hundred milligrams of myometrium was homogenised in 1 mL TRIzol using an IKA Ultra Turrax homogenizer (Cole-Palmer, UK) and extracted using TRIzol plus RNA Maternal request (n) 3 1 3

Medical condition unrelated to pregnancy (n) 3
Poor fetal growth (n) 1 1 1 purification kit (Invitrogen, UK) according to the manufacturer's instructions. The RNA was eluted in 30μL of nuclease-free water and stored at 80˚C. The quantity and purity of RNA was determined using a Qubit RNA BR Assay Kit and fluorometer (Invitrogen, UK) and Nano-Drop 1 ND-1000 spectrophotometer (Thermo Scientific, UK). A 260/ A 280 and A 260 /A 230 absorbance ratios were between 1.8-2.0 and 2.0-2.2 respectively. RNA was treated with TurboDNAse (Invitrogen, UK) to remove any residual genomic DNA. RNA integrity was inspected using an Agilent 2100 Bioanalyser (Agilent Technologies, CA, USA). RNA with an integrity number (RIN) of �7.0 was considered acceptable.

RNA-Seq analysis
cDNA library preparation and sequencing. Total RNA was submitted to the Centre for Genomic Research Facility, University of Liverpool, for RNA-Seq library preparation and sequencing. One microgram of DNA-free total RNA was selected for poly-A using NEBNext Poly (A) mRNA magnetic isolation module (New England Biolabs UK Ltd, UK). RNA-Seq libraries were prepared from the enriched RNA using the NEBNext Ultra Directional RNA library Prep kit for Illumina. Libraries were purified using AMPure XP beads. Each library was quantified using Qubit and the size distribution assessed using the Agilent 2100 Bioanalyser.
Libraries were pooled and assessed by a Qubit dsDNA HS assay kit and qPCR using the Illumina Library Quantification kit from Kapa on a Roche Light Cycler LC48011. Concentrations of template DNA were adjusted to 3nM and denatured with 0.1N NaOH. Denatured cDNA was further diluted to a final concentration of 300pM. Clustering of the DNA templates was performed using HiSeq 3000/4000 paired end cluster Kit with a cBot cluster generation system (Illumina, San Diego, USA) according to the manufacturer's instructions. Libraries were sequenced on an Illumina HiSeq 4000 using sequencing by synthesis chemistry to generate 2 x 150bp paired end reads.
RNA-seq data processing. Sequence libraries for each sample were processed (base called and de-multiplexed) using CASAVA (version 1.8.2, Illumina) to produce sequence data in fastq format. Raw data were filtered to remove Illumina adaptor sequences using Cutadapt version 1.2.1 [37] with option '-O 3' applied to trim adapters from reads if they matched the adapter sequence for 3 bp or more at the 3 0 end. Low quality bases and reads shorter then 20bp were also removed from the data set using Sickle version 1.200.
The trimmed read pairs were aligned using TopHat2 version 2.1.0 [38] against the H. Sapiens genome reference (Ensembl GRCH38, release 93). Default settings were applied, except for the option-g 1. Reads aligning to the reference genome were counted using HTSeq-count version 0.6.1p1 [39].

Differential gene expression analysis
The count data were normalised across libraries to account for differences in library size among samples by using trimmed mean M (TMM) values in edgeR with default parameters. Gene-wise dispersion parameters were estimated and then used for logFC (log 2 fold change) estimating and testing. Differential expression analysis was carried out using edgeR package [40] on all the non-zero count genes identified by TopHat. Differentially expressed genes (DEGs) were extracted by applying the threshold false discovery rate (FDR) of less than 0.05 to adjusted P values, which were generated using Benjamini and Hochberg approach [41].
We performed the following 4 contrasts:

Singleton Term vs Twin Preterm (ST vs. TP)
Contrasts 1 and 2 were performed to examine for differences between singleton and twin pregnancies when matched for gestational age group. Contrasts 3 and 4 were performed to further examine for the effect of stretch on gene expression, with particular emphasis on contrast 3 which was assumed to have the most marked difference in levels of myometrial stretch.

Functional analysis: KEGG pathway and GO term-based analysis
Kyoto Encyclopaedia of Genes and Genomes (KEGG) pathway analysis and Gene Ontology term (GO term) annotation for biological processes, molecular function and cellular components were detected using "gage" and "GOFunction" in R package respectively. P values were adjusted to correct the effects of multiple tests using the Benjamini-Hochberg approach where applicable [41]. Significantly enriched KEGG pathways and GO terms were identified as those with an adjusted P value <0.05.

Data deposition
All sequence data produced in this study have been submitted to ArrayExpress, European Bioinformatics Institute. Raw data files can be accessed under accession number E-MTAB-8233 (https://www.ebi.ac.uk/arrayexpress/experiments/E-MTAB-8233/).

Overview of RNA-Seq data
In total, across the 24 sample cDNA libraries, 1786.7 million reads were sequenced, producing an average of 37.2 million paired-end reads per sample. For all samples, post adapter and quality trimming, the percentage of unpaired reads was less than 0.75%. Unpaired reads were generally shorter than paired reads indicating that they lost more base pairs due to poor quality. Therefore, paired reads post-trimming were better quality than the individual forward or reverse reads. With respect to alignment, most of the reads were successfully aligned (91-93%) using TopHat2 (Table 2). All aligned reads fall in the exon regions of 42,356 genes which represents 72.55% coverage of the current hg38 (version 93) reference transcriptome (58,395 genes in total).

Assessment of variation in count data
To assess clustering of groups, principal component analysis (PCA) of the first two principal components using the log10 count data from all libraries was also performed (Fig 1). This confirmed that sample groups did not cluster into distinct groups and could not be clearly separated. A sample correlation heat map of the correlation of all libraries was also created (Fig 2). This also showed that the correlations are similar both within and between groups. The lowest correlation amongst groups was >0.92. Together this indicted that the number of differentially expressed genes detected would be small.

Identification of differentially expressed genes
When matched for gestational age (contrasts 1 and 2), we found only 3 genes which were differentially expressed between singleton and twin myometrium (Table 3). These genes correspond to: ADRA1D, alpha 1D adrenergic receptor, and DPY19L2, dpy-19-like 2 which were up-regulated in preterm twin pregnancy compared to preterm singleton pregnancy and HLA-DQA2, major histocompatibility complex class II, DQ Alpha 2, which was down-regulated in twin pregnancy at term compared to term singleton pregnancy.
In both gestation groups (preterm and term), the combined birthweight of the neonates in twin pregnancies was significantly greater than the neonatal birthweight in singleton pregnancies ( Table 1). The mean birthweight for the twin pregnancy group at term however, was more than double that of the mean birthweight of the singleton preterm group (5525g and 2726g respectively, Table 1). The additional stress placed on the myometrium in twin pregnancy resulting from the increased intrauterine volume has been suggested to contribute to higher incidence of preterm birth in this group. Using birthweight as an estimate of the level of uterine wall stress and hence a surrogate for investigating the effects of crude stretch on myometrial gene expression, we next compared the transcriptome data for twin pregnancies at term (TT) with singleton preterm pregnancies (SP) where the differences in birthweight and hence volume would be at its greatest.
In total, the expression of 99 transcribed elements, including protein-coding transcripts, pseudogenes and long intergenic noncoding RNAs (lincRNAs) were significantly different in twin myometrium compared with singleton myometrium. Of these, 24 were at higher levels in the twin myometrium and 75 were at lower levels in the twin myometrium. Of the 24 transcripts expressed at a higher level in twin myometrium, 14 were known protein-coding genes. The remaining 10 genes contained 4 pseudogenes, 3 lincRNAs and 3 antisense genes. Of the 75 downregulated genes in twins, 70 were known protein coding genes. The remaining 5 gene biotypes included 4 antisense genes and one IG gene. Complete lists of up-and down-  Tables 4 and 5 respectively. The entire list of significantly DEGs is presented in S1 Table. Cluster analysis of the DEGs for all contrasts is shown in Fig 3. The analysis identified differences in expression of a number of genes associated with inflammation and immune response, including: HLA-DQA2, IL19, IL2RB, IL21R, CXCR6, CCL5, XCL2 and LBP as well as CD7, GATA3, TNFRS18, CDE3, GNLY and GZMA which are   genes involved in T cell maturation, differentiation and regulation of interactions with target cells. These were all identified as being down-regulated (log 2 fold change �-1, p<0.05) in twin term myometrium. Other down-regulated genes included 3β-hydroxysteroid dehydrogenase (HSD3B1) and steroidogenic acute regulatory protein (STAR) which are involved in steroid hormone biosynthesis, including the catabolism of cholesterol to pregnenolone and the conversion of pregnenolone to progesterone. Genes associated with regulation of contraction such as KCNG1 and KCNH1 members of the voltage-sensitive potassium (K v ) channel family, the interstitial cells of Cajal-like cell marker, c-kit (CD177) and phosphodiesterase 11A (PDE11A) which is involved in the downregulation of cAMP and cGMP signalling were also identified as being differentially expressed in twin term myometrium.
Comparing term singleton myometrium with preterm twin myometrium (contrast 4) however, only identified one gene which was differentially expressed. This corresponded to, carbonic anhydrase 1 (CA1) which was found to be up regulated in preterm twin myometrium compared to term singleton myometrium.

GO functional enrichment and KEGG pathway analysis
To systematically determine functional analyses and canonical pathways that the DEGs might be involved in, we planned to perform GO-term and KEGG pathway enrichment on those genes identified as being significantly up-or down-regulated in twin myometrium. As the number of DEGs identified in contrasts 1, 2 and 4 was small, the DEG data generated was not suitable for use in pathway analysis. We therefore focussed on contrast 3: Singleton preterm versus twin term myometrium.
Genes were assessed for biological process, molecular function and cellular component GO-term enrichment and KEGG pathway enrichment using "GOFunction" and "gage" in R software. The genes identified were assigned to 24, 8 and 16 functional groups in the different GO categories respectively (Tables 6-8). The top GO-term biological processes identified as being altered in twin term myometrium compared to preterm singleton myometrium included 'immune response,' 'regulation of immune response' and 'inflammatory response' (Table 6). Whilst molecular function and cellular component ontologies were associated mostly with 'receptor binding' and the 'plasma membrane' (Tables 7 and 8).
Active KEGG pathways also included cytokine-cytokine receptor interaction, chemokine signalling, toll like receptor signalling and nod-like receptor signalling which are all important inflammatory pathways mapped to the 'immune system' and 'signalling molecules and interaction' biological categories (Table 9).

Discussion
In this study, we examined for the first time, the human myometrial transcriptome from singleton and twin pregnancies at two gestational time points (preterm and term) using high throughput RNA sequencing. When matched for gestational age, we find only a few genes which are differentially expressed between groups which implies that there is little or no difference in the myometrial transcriptome of singleton and twin pregnancies, at least when length of pregnancy is controlled for. Our unbiased statistical analysis of the sequenced transcriptomes of the 24 analysed samples in the form of principal component analysis was not able to clearly distinguish between the 4 clinical subgroups; singleton preterm, twin preterm, singleton term and twin term, further reflecting that there is little difference between the transcriptomes of these samples. Previous studies employing RNA-Seq in human myometrium, have successfully examined changes in genes associated with labour onset in singleton myometrium in which there is a clear difference in the transcriptome once labour has been initiated [34,36]. The DEGs identified in these RNA-Seq studies also broadly reflect those of earlier microarray studies. RNA-Seq also has the advantage of increased sensitivity and is not limited to the complement of preselected probes present on an array. Hence, we are confident that a high throughput RNA sequencing approach was suitable for use in this study.
Previous data relating to changes in gene and protein expression in twin pregnancy has been restricted to studies of a single or few genes [10,22,23] particularly focussed around myometrial stretch. Whilst this is useful for pinpointing the role of specific genes and proteins in a particular pathway, it does not give an overall unbiased view or qualitative assessment of the global changes in mRNA expression within or between tissues at a given time, which we aimed to identify here. During human pregnancy, the uterus grows and expands to accommodate the developing fetus. Early in gestation this is mainly via hyperplasia (increased number of uterine cells) and a small contribution from increased cell size (hypertrophy). Towards late pregnancy, growth of the fetus imposes tension on the uterine wall, inducing a number of biochemical and molecular changes in the myometrial cells [42,43]. This stretch is thought to be the major stimulus for uterine growth [44] and is a continual process allowing for uterine remodelling. Stretch pressures placed on the uterine wall by the growing fetuses in animal models of pregnancy and labour has been shown to increase the expression of myometrial contraction-associated protein (CAP) genes including connexin 43 (gap junction protein) and the oxytocin receptor (OTR) as well as to changes in the number of extracellular matrix proteins and activation of focal adhesions [45]. Application of mechanical stretch to human myometrial cells in vitro is also associated with up-regulation of OTR, COX2 and secretion of cytokines and IL-8 [11][12][13], further supporting a role for myometrial stretch contributing to myometrial activation.
It is not unreasonable to expect changes in myometrial tension in twin pregnancies given they tend to have a 60% increase in overall fetal birthweight. The hypothesis surrounding myometrial stretch and PTL onset is that PTL may reflect an inability of the uterus to adapt to this continuous increase in volume. In multiple pregnancies therefore (and those complicated by polyhydramnios), whilst stretch is a normal process and an important stimulus for uterine growth, stretch may be excessive. Hence, twin pregnancies which are 'stretched' to a greater extent at an equivalent time point in gestation compared to singleton pregnancies, are associated with an increased risk of preterm labour and delivery [46]. Undoubtedly, duration of pregnancy is shorter with increasing numbers of fetuses [47]. A recent large Swedish retrospective cohort study also supports uterine distention playing a causal role in determining the Table 6. Significant Gene Ontology terms identified for biological processes in twin term myometrium compared to singleton preterm myometrium.

GO-term name: Biological Process
Number of genes linked to this term % of total genes in the list ( timing of birth in uncomplicated spontaneous deliveries, as well as a role for increased uterine load, maternal height and uterine size [48]. In examining tissues which were matched for gestational age and for other characteristics (as reported here), one may therefore have expected to see changes in gene expression resulting from the assumed increase in pressure on the uterine wall in twin pregnancy.
Our data showing no change in term samples is in agreement with a study by Lyall et al., [23] who investigated the expression of key regulators of myometrial function including, connexin 43 and 26, prostanoid receptors, EP1-EP4 and Gαs, in term (~38 weeks) singleton and twin myometrium and found no difference in expression between the two groups. They concluded that the mechanisms whereby stretch can activate myometrial contraction are complex and hence other factors must be involved. Our gestational age-matched cohort RNA-Seq data however, did not reveal any other candidate genes that warrant further investigation: The ADRA1D gene encodes alpha-1D-adrenergic receptor which is a member of the G proteincoupled receptor superfamily. Adrenergic receptors activate mitogenic responses and regulate growth and proliferation of many cells hence its role in myometrial cell growth in twins could be supported. HLA-DQA2 is involved in antigen processing and presentation of exogenous peptide antigens for CD4 T-cell recognition via MHC classs II molecules which are expressed in antigen presenting cells (B lymphocytes, dendritic cells, macrophages). Hence its expression here probably arises from other cells i.e. immune cells, present within the myometrial biopsy.
To further explore the potential effect of stretch on the myometrial transcriptome, we next examined for differences in gene expression between singleton preterm and twin term samples as in terms of birthweight, this reflects the greatest difference between groups (more than double) and hence the greatest assumed change in uterine wall stress and myometrial stretch. In doing so, we detected 99 DEGs of which 74 were known protein coding genes and included changes in pro-inflammatory genes, genes associated with regulation of myometrial contraction and genes involved in the de-novo synthesis of steroids such as pregnenolone and progesterone which are important for the maintenance of pregnancy. Interestingly these were all down-regulated in the twin group.
We acknowledge however, that the exact increases in uterine wall tension cannot be determined and note that whether the significant increase in intrauterine volume in twin pregnancies also affects uterine wall tension is questioned [49]. Is it therefore difficult to say whether these genes reflect the crude/long term effects of stretch on myometrial gene expression. The changes identified could also be due to changes in gestational age between the two groups. Whilst it was not the focus of our study, we addressed this by also performing a subset analysis to identify DEGs between preterm and term groups (data not shown). We found 32 DEGs between the singleton preterm and term groups, but, only 12 of these (~10%) overlapped with the DEGs from contrast 3. We found no changes between twin preterm and term samples. Hence, we do not think that differences in gestational age contribute extensively to the DEGs identified between preterm singleton and term twin myometrium. Our analysis of Gene Ontology terms relating to biological processes and KEGG pathways between singleton preterm and twin term samples also revealed alterations in immune response, regulation of immune response and inflammation, as well as cytokine and chemokine signalling. The inflammatory genes identified however, were down-regulated in our twin group which is contradictory to stretch inducing pro-inflammatory gene expression as previously reported from in vitro studies examining the effects of stretch on human myometrial smooth muscle cells [10,13,14,50]. One explanation for this could be that in vivo, the stretch applied is not acute and the myometrium may be responding to and compensating for the increasing stretch over time. The effects of chronic (long term) stretch such as anticipated in twins here is therefore likely to be different to the effects of the more acute mechanical stretch as used in studies investigating the effects of stretch in vitro.
Changes in protein function with increased stretch rather than levels of expression may also be playing a role and hence would not be detected in our study. For example, Aye et al, recently suggested that stretch of myometrium induces a change in the conformation of the transmembrane domains of the oxytocin receptor which facilitates signal transduction and leads to enhanced contractility even in the absence of oxytocin [51]. Hence in multiple pregnancy, stretch-induced activation of the OTR rather than increased OTR expression, may contribute to the higher rates of preterm birth in these women.

Study limitations
One of the main limitations of this study is that samples were obtained from women undergoing 'elective' pre-labour CS. We chose to examine samples from these women with the intention of being able to study changes in gene expression between groups independently of changes that may occur with labour onset. We felt that using samples collected at the time of labour (preterm or term), irrespective of the initial mechanistic trigger, could be too biologically similar to distinguish changes in expression between sample groups. For example, the activation of similar final common pathways leading to labour such as inflammation, may mask any significant detectable changes between groups. To this end we also chose samples from women with similar base characteristics including gestational age (within gestation group), parity, height, weight and BMI and maternal age to further reduce non labour-associated variation between samples. In addition, therapeutically, any potential gene marker of PTL would need to be identifiable prior to the onset of labour to allow the greatest opportunity to intervene to prevent preterm delivery.
In choosing women pre-labour however, we were unable to fully match samples by indication for CS. This is because elective preterm deliveries are usually complex and performed to protect the health of the mother and/or fetus from continuing pregnancy, whilst elective CS at term are most likely due to fetal malposition, previous CS or complications in a prior delivery. In this study, we minimised the impact of indication for CS on study outcomes by only including preterm women with medical issues if their condition was unrelated to pregnancy, and excluded any woman on medication e.g. for diabetes or hypertension (all gestations).
We also hypothesised that this study would enable us to isolate changes in the transcriptome resulting from the extra biomechanical stress or 'stretch' placed on the myometrium, as anticipated in twin pregnancy. Whilst little or no changes were observed in our groups when matched for gestational age, our analysis of singleton preterm myometrium verses twin term myometrium (where the difference in neonatal birthweight is most significant) did detect some DEGs. However, we cannot rule out the effects of other factors e.g. hormonal differences (fetal or maternal) in multiple gestations being responsible for the differences observed. As our samples were from pre-labour women, we also do not know how 'far off' from labour the patients were and therefore cannot control for changes in gene expression which may be associated with the myometrium transitioning towards labour onset.
It is disappointing that the results of this study do not add much to our understanding of the different mechanistic triggers for preterm (and term) labour in twins verses singletons, should one exist. The high correlation between our samples (>0.92) and lack of sample group separation by principle component analysis supports that the real number of DEGs between groups is small. Increasing study sample size is therefore unlikely to increase the identification of DEGs between our gestational age-matched cohorts. Application of more stringent criteria for gene filtering (i.e. removing genes with zero counts in >25% of samples, did not significantly alter DEG number (data not shown) nor did it provide further sensible functional interpretation of our findings.
Examining expression from in-labour samples may yield greater differences. Hence, future studies should aim to examine the myometrial transcriptome changes associated with preterm labour onset in both singletons and twins, comparing samples from those in not in labour to those where spontaneous labour has begun. For example, Johnson's group comparing twins in preterm labour versus no labour showed that protein levels of inflammatory cytokines and chemokines including TNFα, Il6, IL8 and CCL2 were increased in myometrial tissues [10] but they did not examine singleton tissues. By their nature, however, preterm in-labour tissues are more difficult to obtain since most women delivering preterm will do so vaginally. In addition, further examination of changes in gene expression according to parity (i.e. comparing nulliparous versus parous women) would also be informative, since parity is known to affect pregnancy and labour-specific processes.
In conclusion, twin pregnancy is associated with higher rates of preterm labour, possibly attributed to increased myometrial stretch and uterine over-distention. We find little or no difference in the transcriptome of singleton and twin pregnancy myometrium when matched for gestational age. In this regard we are therefore unable to conclusively add to previously suggested mechanisms, including uterine stretch, which are responsible for premature activation of myometrium in multiple pregnancy. Factors derived from other tissues such as the placenta, amnion and the fetus itself may be important, as well as examining samples from women undergoing spontaneous preterm labour.