Objective
The purpose of this study was to determine whether cervical shortening of a ripe cervix at term is associated with changes in the cervical transcriptome.
Study Design
Sonographically measured cervical lengths and biopsy specimens were obtained from 19 women at term who were not in labor with a ripe cervix. Affymetrix HG-U133 Plus 2.0 arrays (Affymetrix Inc, Santa Clara, CA) were used. Gene expression was analyzed as a function of cervical length. Gene Ontology, pathway analyses, quantitative real-time reverse transcription–polymerase chain reaction, and immunohistochemistry were performed.
Results
Cervical length shortening was associated with differential expression of 687 genes. Fifty-four biologic processes, 22 molecular functions, and 9 pathways were enriched. Quantitative real-time reverse transcription–polymerase chain reaction analysis confirmed differential expression of 13 genes. Bone morphogenetic protein-7, claudin-1, integrin beta-6, and endometrial progesterone–induced protein messenger RNA, and protein expressions were down-regulated with cervical shortening.
Conclusion
Sonographic cervical shortening in patients at term who are not in labor with a ripe cervix is associated with changes in the uterine cervix transcriptome. The epithelial-mesenchymal transition may participate in the mechanism of cervical shortening at term.
The uterine cervix plays a central role in the maintenance of normal pregnancy and in parturition. Before the onset of labor (both at term and preterm) changes in the uterine cervix have been reported with the use of digital examination and sonography. A sonographic short cervix has been used to predict the onset of labor and successful induction of labor and is considered to be the most powerful predictor of spontaneous preterm birth. Shortening of the cervix that is detected by transvaginal sonography has been used to identify patients who may benefit from cerclage placement or progesterone treatment. Despite its importance, little is known about the molecular mechanisms that are responsible for cervical shortening.
Most studies about the physiology, biophysical properties, biochemistry, and cellular composition of the uterine cervix have been conducted by hypothesis-driven research. Observations that have been derived from these studies have contributed to the understanding of the mechanisms that are responsible for cervical ripening and remodeling. High-dimensional biology allows an unbiased description of the changes that occur with cervical shortening without the use of biologic assumptions. The aim of this study was to use transcriptomics to identify changes in gene expression as a function of cervical shortening.
Materials and Methods
A cross-sectional study was performed in women at term and not in labor with a ripe cervix who were scheduled to have an elective cesarean section. A ripe cervix was defined as a Bishop score of >5. Women were asked to participate in a protocol in which transvaginal ultrasound examination of the cervix was performed before delivery; patients were consented to have a cervical biopsy while under anesthesia during cesarean section. Inclusion criteria were (1) term gestation (≥37 weeks), (2) no prostaglandin or oxytocin administration, (3) negative Neisseria gonorrhoeae and Chlamydia trachomatis determined by examination of cervical secretions, and (4) a normal Papanicolaou smear test. Patients were excluded if they had evidence of histologic chorioamnionitis. Eight patients from a previous study were included in the microarray analysis, and 11 additional patients were added for the quantitative real-time reverse transcription–polymerase chain reaction (qRT-PCR) analysis. This study was approved by the institutional review board of Wayne State University, and all patients provided written, informed consent. Half centimeter biopsy specimens were obtained transvaginally from the anterior lip of the cervix at the 12 o’clock position and immediately snap frozen in liquid nitrogen or placed in RNAlater (Ambion Inc, Austin, TX) and stored at –70°C. This procedure has been used extensively by investigators in the United States, Europe, and other continents. The demographic and clinical data, obstetric and gynecologic history, and pregnancy outcome were extracted from medical records.
Microarray analysis
Total RNA was extracted, and the transcriptome was profiled with Affymetrix HG-U133PLUS 2.0 arrays (Affymetrix Inc, Santa Clara, CA) as previously described. Gene expression data preprocessing that included background correction and normalization was performed with the robust multichip analysis algorithm implemented in the Affy package of Bioconductor (Affymetrix Inc). Further analysis was performed on the probesets that were determined to be “present” (4/8 samples) with the Affymetrix MAS5.0 detection call method (Affymetrix Inc). A linear model was fit on the log-transformed gene expression levels as a function of the cervical length of all patients; adjustment was made for maternal age and number of previous vaginal deliveries. The estimated coefficient for the cervical length (slope) for each probe represents an overall (average) rate of increase or decrease of the gene expression as a function of cervical length. The use of cervical length as a continuous variable rather than a dichotomous variable was preferred in this analysis because it provides greater power and alleviates the ambiguity of choosing a cutoff of cervical length.
The limma software package that implements a moderated t -statistic was used to derive a probability value for each probeset. A false discovery rate (FDR) of <0.2 and fold change of >1.5 per 1 cm of change in cervical length was used to select differentially expressed genes.
To gain biologic insight from the resulting differentially expressed genes, a Gene Ontology analysis was conducted. The Gene Ontology analysis provides a hierarchically structured classification of genes under the broad categories “molecular function,” “biological process,” and “cellular component.” The Gene Ontology enrichment analysis was conducted with the GOstats library of bioconductor.
Pathway Analysis was performed on the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway database with the use of both an enrichment analysis and the Signaling Pathway Impact Analysis. An FDR of <0.05 was used for Gene Ontology and Pathway Analysis. For the gene expression analysis, an FDR threshold of 0.2 was used because it was combined with an additional layer of stringency based on the magnitude of change (at least 1.5-fold per cm). Thus, among all probes with the corrected FDR of <0.2, only those that changed at least 1.5-fold per cm of cervical length were selected as differentially expressed. Therefore, the overall FDR is expected to be <0.2. For the Gene Ontology and Pathway Analyses, only the FDR was used as the criteria to select significant Gene Ontology terms and pathways; therefore, a more stringent threshold of 0.05 was used.
qRT-PCR analysis
Selected significant genes were tested by qRT-PCR to confirm the microarray results. The Biomark System (Fluidigm, San Francisco, CA) was used to perform high-throughput qRT-PCR confirmation. For this system, specific target amplification of complementary DNA was performed. Briefly, 0.2X pool of specific gene expression assays (Applied Biosystems, Foster City, CA) was used as the source of primers. Preamplification reactions contained 1.25 μl complementary DNA, 2.5 μl TaqMan PreAmp master mix (Applied Biosystems), and 1.25 μl pooled assay mix. The reaction was performed with a thermal cycler for 14 cycles at 95°C for 15 seconds and 60°C for 4 minutes. After the cycling, the reaction was diluted 1:5 by double distilled water to a final volume of 25 μl. Fluidigm 96.96 Dynamic Array chip was used to perform the next step qRT-PCR assays. The 96.96-array chip was primed in an integrated fluidic circuit controller. After the priming, 2.5 μl 20X TaqMan gene expression assays (Applied Biosystems) were mixed with 2.5 μl 2X assay loading reagent (Fluidigm) and loaded into the assay inlet on the 96.96-array chip; 2.25 μl preamplified complementary DNA was mixed with 2.5 μl TaqMan Universal PCR master mix (Applied Biosystems) and 0.25 μl 20X GE sample loading reagent (Fluidigm) and loaded into the sample inlet on the chip. The chip was returned to the integrated fluidic circuit controller for loading. After the samples and assays were loaded, the chip was placed into the Biomark System to run the reactions. The cycle threshold (Ct) value of each reaction on the chip was obtained with the Fluidigm RT-PCR analysis software. The surrogates for log 2 gene expression levels, –ΔCt values, were obtained from the target and reference gene Ct values. As with the microarray data, a linear model was used to further test the association between the expression level and the cervical length, while adjustment was made for maternal age and the number of previous vaginal deliveries.
Immunohistochemistry
Immunohistochemistry for endometrial progesterone–induced protein (EPIP), claudin-1, and bone morphogenetic protein-7 (BMP-7) was performed on 5-μm tissue sections of frozen human cervical tissue that were cut and placed on positively charged slides. Primary antibodies included rabbit polyclonal anti-human EPIP (clone 133929; 1:50 dilution; Santa Cruz Biotechnology Inc, Santa Cruz, CA), mouse monoclonal anti-human claudin-1 (clone ab56417; 1:400 dilution; Abcam Inc, Cambridge, MA), and mouse monoclonal anti-human BMP-7 (clone 164311; 1:40 dilution; R&D Systems, Minneapolis, MN) antibodies. The primary antibodies were applied and incubated for 1 hour, which was followed by washing and incubation with the biotinylated secondary antibody for 30 minutes at room temperature. Detection was performed with AEC counterstaining with Mayer hematoxylin followed by mounting.
Results
Demographic and clinical characteristics of the study population are depicted in Table 1 . The median cervical length of the 19 women before cesarean delivery was 23 mm (range, 5–49). Approximately one-half of these patients (52%) had a sonographic cervical length of ≤25 mm. Of those, 2 patients had a cervical length of <1.0 cm, and 7 patients had a cervical length of <2.0 cm. Of the 19 patients, 74% of the women (14/19) were African American; 95% of the women (18/19) did not smoke, and 84% of the women (16/19) had “elective repeat” as the indication for their cesarean section delivery.
| Variable | Subjects included in the microarray analysis (n = 8) a | Subjects included in the quantitative real-time reverse transcription–polymerase chain reaction analysis (n = 19) a |
|---|---|---|
| Age, y | 28 (20–40) | 28 (20–40) |
| Parity, n | 2 (0–5) | 2 (0–5) |
| Previous vaginal deliveries, n | 0 (0–4) | 0 (0–4) |
| Cervical length, mm | 27 (16–49) | 23 (5–49) |
| Bishop score | 6.5 (5–9) | 6 (5–8) |
| Cervical dilation, cm | 1 (0–3.5) | 1 (0–3.5) |
| Gestational age at delivery, wk | 39 (37–40) | 39 (37–40) |
Microarray analysis
Eight hundred thirty five probesets corresponding to 687 genes were differentially expressed as a function of cervical length after controlling for maternal age and number of prior vaginal deliveries. The expression of 226 genes was down-regulated and of 461 genes was up-regulated as cervical length shortened. The top 50 genes that were differentially expressed are listed in Table 2 . An additional analysis was conducted to ensure that there was no confounding between the body mass index and cervical length. Inclusion of the body mass index as a covariate in the model did not substantially affect the estimated gene expression changes (direction and amount of change and gene rank) as a function of cervical length.
| Entrez gene | Symbol | Gene name | Fold change with each 30-mm change in cervical length | Fold change with each 10-mm change in cervical length | Direction of change with cervical shortening |
|---|---|---|---|---|---|
| 1828 | DSG1 | Desmoglein 1 | 346.36 | 7.02 | ↓ |
| 4604 | MYBPC1 | Myosin binding protein C, slow type | 264.44 | 6.42 | ↑ |
| 3294 | HSD17B2 | Hydroxysteroid (17-beta) dehydrogenase 2 | 187.11 | 5.72 | ↑ |
| 3240 | HP | Haptoglobin | 171.78 | 5.56 | ↑ |
| 5172 | SLC26A4 | Solute carrier family 26, member 4 | 115.07 | 4.86 | ↑ |
| 89778 | SERPINB11 | Serpin peptidase inhibitor, clade B (ovalbumin), member 11 (gene/pseudogene) | 113.50 | 4.84 | ↓ |
| 6374 | CXCL5 | Chemokine (C-X-C motif) ligand 5 | 79.76 | 4.30 | ↑ |
| 7134 | TNNC1 | Troponin C type 1 (slow) | 76.59 | 4.25 | ↑ |
| 276 | AMY1A | Amylase, alpha 1A (salivary) | 69.49 | 4.11 | ↑ |
| 1080 | CFTR | Cystic fibrosis transmembrane conductance regulator (ATP-binding cassette sub-family C, member 7) | 63.97 | 4.00 | ↑ |
| 28823 | IGLV1-44 | Immunoglobulin lambda variable 1-44 | 63.16 | 3.98 | ↑ |
| 8424 | BBOX1 | Butyrobetaine (gamma), 2-oxoglutarate dioxygenase (gamma-butyrobetaine hydroxylase) 1 | 61.50 | 3.95 | ↓ |
| 1365 | CLDN3 | Claudin 3 | 59.12 | 3.90 | ↑ |
| 27324 | TOX3 | TOX high mobility group box family member 3 | 57.74 | 3.87 | ↑ |
| 27324 | TOX3 | TOX high mobility group box family member 3 | 57.19 | 3.85 | ↑ |
| 155465 | AGR3 | Anterior gradient homolog 3 (Xenopus laevis) | 56.45 | 3.84 | ↑ |
| 3500 | IGHG1 | Immunoglobulin heavy constant gamma 1 (G1m marker) | 55.31 | 3.81 | ↑ |
| 27324 | TOX3 | TOX high mobility group box family member 3 | 53.46 | 3.77 | ↑ |
| 3535 | IGL@ | Immunoglobulin lambda locus | 52.66 | 3.75 | ↑ |
| 9185 | REPS2 | RALBP1 associated Eps domain containing 2 | 50.75 | 3.70 | ↑ |
| 115111 | SLC26A7 | Solute carrier family 26, member 7 | 49.00 | 3.66 | ↑ |
| 79949 | C10orf81 | Chromosome 10 open reading frame 81 | 48.18 | 3.64 | ↑ |
| 148808 | MFSD4 | Major facilitator superfamily domain containing 4 | 44.85 | 3.55 | ↑ |
| 64699 | TMPRSS3 | Transmembrane protease, serine 3 | 44.69 | 3.55 | ↑ |
| 5169 | ENPP3 | Ectonucleotide pyrophosphatase/phosphodiesterase 3 | 43.43 | 3.52 | ↑ |
| 121838 | LOC121838 | Hypothetical LOC121838 | 42.39 | 3.49 | ↑ |
| 9185 | REPS2 | RALBP1 associated Eps domain containing 2 | 40.68 | 3.44 | ↑ |
| 7103 | TSPAN8 | Tetraspanin 8 | 40.40 | 3.43 | ↑ |
| 9173 | IL1RL1 | Interleukin 1 receptor-like 1 | 39.76 | 3.41 | ↑ |
| 3170 | FOXA2 | Forkhead box A2 | 39.38 | 3.40 | ↑ |
| 400169 | DKFZp451A211 | DKFZp451A211 protein | 38.79 | 3.39 | ↑ |
| 10551 | AGR2 | Anterior gradient homolog 2 (Xenopus laevis) | 38.72 | 3.38 | ↑ |
| 79679 | VTCN1 | V-set domain containing T cell activation inhibitor 1 | 37.25 | 3.34 | ↑ |
| 643187 | LOC643187 | Hypothetical LOC643187 | 36.96 | 3.33 | ↓ |
| 2952 | GSTT1 | Glutathione S-transferase theta 1 | 36.40 | 3.31 | ↑ |
| 51237 | MGC29506 | Hypothetical protein MGC29506 | 35.84 | 3.30 | ↑ |
| 3851 | KRT4 | Keratin 4 | 35.29 | 3.28 | ↓ |
| 387914 | SHISA2 | Shisa homolog 2 (Xenopus laevis) | 34.95 | 3.27 | ↑ |
| 3535 | IGL@ | Immunoglobulin lambda locus | 34.75 | 3.26 | ↑ |
| 5803 | PTPRZ1 | Protein tyrosine phosphatase, receptor-type, Z polypeptide 1 | 34.68 | 3.26 | ↓ |
| 57864 | SLC46A2 | Solute carrier family 46, member 2 | 34.14 | 3.24 | ↑ |
| 3535 | IGL@ | Immunoglobulin lambda locus | 33.79 | 3.23 | ↑ |
| 54916 | C14orf101 | Chromosome 14 open reading frame 101 | 33.06 | 3.21 | ↑ |
| 23231 | KIAA0746 | KIAA0746 protein | 32.41 | 3.19 | ↑ |
| 155006 | TMEM213 | Transmembrane protein 213 | 32.23 | 3.18 | ↑ |
| 1010 | CDH12 | Cadherin 12, type 2 (N-cadherin 2) | 31.53 | 3.16 | ↑ |
| 92196 | DAPL1 | Death associated protein-like 1 | 31.06 | 3.14 | ↓ |
| 64067 | NPAS3 | Neuronal PAS domain protein 3 | 30.18 | 3.11 | ↑ |
| 608 | TNFRSF17 | Tumor necrosis factor receptor superfamily, member 17 | 30.08 | 3.11 | ↑ |
| 199920 | C1orf168 | Chromosome 1 open reading frame 168 | 29.97 | 3.11 | ↑ |
Stay updated, free articles. Join our Telegram channel
Full access? Get Clinical Tree
