Materials and Methods
This was a case-control genetic association study. Women of non-Hispanic European ancestry with at least 1 previous documented singleton, nonanomalous SPTB who delivered at <34.0 weeks’ gestation who received 17P in at least 1 subsequent pregnancy were recruited prospectively from a consultative Preterm Birth Prevention Clinic from 2008-2010 at Intermountain Medical Center (Salt Lake City, UT). The treatment of women in the prematurity prevention clinic has been described previously. All women provided written informed consent, and this study was approved by the institutional review boards at Intermountain Healthcare and the University of Utah.
Women with a history of either idiopathic spontaneous onset of contractions and cervical dilation or preterm premature rupture of membranes (PPROM) in a previous pregnancy were considered to have a previous SPTB and were included. We excluded women with a known or suspected cause of SPTB, which included women who experienced SPTB after polyhydramnios, within 2 weeks of amniocentesis, because of hypertensive disorders that included preeclampsia or were related to abdominal trauma. Women with known uterine anomalies or a history of treatment for cervical dysplasia with cryotherapy, loop electrosurgical excision procedure, or cervical conization were also excluded.
Women were classified as a 17P responder or nonresponder based on response to 17P. Specifically, the difference between the earliest delivery gestational age because of SPTB without 17P and the delivery gestational age with 17P was calculated and was termed the “17P effect.” If a woman had multiple pregnancies that had been treated with 17P after her initial SPTB, the 17P effects from each individual pregnancy were averaged to generate an overall 17P effect. Women with an overall 17P effect of ≥3 weeks (ie, the individual’s pregnancy or pregnancies that had been treated with 17P delivered at least 3 weeks later compared with the gestational age of the earliest SPTB without 17P treatment) were considered 17P responders. Women with a negative overall 17P effect and those with an overall 17P effect of <3 weeks were classified as nonresponders. Demographic data were compared between responders and nonresponders with the Student t test and Fisher exact test, as appropriate (version 12.1; StataCorp LP, College Station, TX).
DNA was extracted from stored buffy coats. All extracted DNA underwent quality control with a spectrophotometer (NanoDrop Products, Wilmington, DE) reading and evaluation on a 1% agarose gel before genomic library construction. Genomic libraries were then constructed, and samples underwent additional quality control measures that included quantitative polymerase chain reaction quantitation of library concentration with primers (Illumina Inc, San Diego, CA) and evaluation of the library on an Agilent Bioanalyzer DNA 1000 chip (Agilent Technologies Inc, Santa Clara, CA). A PhiX control library (Illumina Inc) was spiked into each lane at a concentration that represented approximately 0.5% of the reads. This platform targeted approximately 180,000 protein-coding exons, in approximately 20,000 genes, for capture. Whole exome sequencing was then performed at The University of Utah Huntsman Cancer Institute’s Microarray Core Facility with Illumina HiSeq2000 (Illumina Inc) technology. We indexed 4 samples per lane, with a goal of approximately ×40-50 average depth of coverage per sample.
Sequences were then called simultaneously on all samples with the University of Utah Department of Human Genetics variant-calling software pipeline. Paired-end 101 base pair fastq reads were aligned to the reference genome (b37) with the Burrows-Wheeler aligner software. Additional processing that included sorting, mate-fixing, and duplicate read removal was performed with Samtools and Picard Tools. Insertion and deletion realignment and base recalibration was performed with the Genome Analysis Tool Kit (Broad Institute, Cambridge, MA). Processed call-ready BAM files were called jointly with the Unified Genotyper (Broad Institute). Raw genotypes were evaluated and filtered with the Variant Quality Score Recalibrator that is provided in the Genome Analysis Tool Kit package.
Individual exomes were analyzed for evidence of population stratification with Eigenstrat software (Harvard School of Public Health, Boston, MA). Single nucleotide polymorphisms with a minor allele frequency of <0.05 and/or a strong deviation from Hardy-Weinberg equilibrium ( P < .00001) were removed. Single nucleotide polymorphisms were also filtered to remove all single nucleotide polymorphisms with pairwise linkage disequilibrium of r 2 >0.2. Population stratification was then analyzed among the remaining subjects, and individual outliers were excluded from further analysis.
The remaining exomes were compared between 17P responders and nonresponders using the Variant Annotation, Analysis, and Search Tool (VAAST). VAAST is a publicly available probabilistic search tool for the identification of disease-causing variants. VAAST scores coding variants that are based on the allele and amino acid substitution frequencies differences between case and control genomes have been demonstrated to be effective in the identification of causative disease alleles both in cases of rare variants in rare disease and combinations of rare and common variants in common disease. The VAAST analysis produced a “raw” list of genes that was prioritized by the likelihood of allelic differences between 17P responders and nonresponders. We reported raw probability values (uncorrected for multiple comparisons) for our VAAST analysis gene list.
Next, the genes that were obtained from the VAAST analysis were compared with those potential candidate genes. Using the online Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway database, we compiled a list of possible candidate genes that we suspected to be involved in either the pathogenesis of prematurity (ie, genes in pathways that affect inflammatory response, myometrial contraction/relaxation, oxidative stress, coagulation, and complement, calcium signaling) or the metabolism of 17P (ie, genes in pathways that involve steroid receptors, drug metabolism, steroid degradation). This list of candidate genes was compared with the raw list of VAAST genes.
In the next step of our pathway analysis, we again used the raw VAAST list of genes and selected the top 2.5% of genes with the highest VAAST scores from that list. Each gene on this raw list was classified by the online Protein Analysis THrough Evolutionary Relationships (PANTHER) system into known gene ontology molecular functions and biologic processes. The percentage of genes within each molecular function/biologic process category was compared with an online referent population to search for areas of over- and under-representation with the use of the binomial test that is available through PANTHER tools online. A Bonferroni-corrected probability value < .05 for the binomial test was considered significant.
Results
Fifty-six women met initial inclusion criteria. All were of self-reported European ancestry. Of these, sequencing analysis failed in 2 women (average sequencing depth of coverage <1). On initial Eigenstrat analysis, population differentiation between cases and control subjects was not significant ( P > .09). However, 4 samples deviated substantially from the main cluster of points. Removal of these 4 samples produced a much less stratified data set with little genome-wide differentiation between cases and control subjects ( P > .35). Thus, the final cohort consisted of 50 women (41 responders and 9 nonresponders). All remaining samples met genotype quality filters. The average depth of exome coverage was 51 ± 18 base pairs (range, 13.4–102.6 base pairs).
Demographic and previous pregnancy characteristics were similar between responders and nonresponders with regard to parity, number of preterm births before the studied gestation, cervical insufficiency history, and PPROM history ( Table 1 ). Four women (3 responders, 1 nonresponder) had a history of a cervical laceration that was related to delivery; this was not significantly different between groups ( P = .56). This group of women was generally healthy. None of the women had a history of type I or II diabetes mellitus or chronic renal disease. Four women (8%) had a history of chronic hypertension; none of the women experienced preeclampsia or required preterm delivery because of worsening hypertension.
| Variable | 17P responders (n = 41) | 17P nonresponders (n = 9) | P value |
|---|---|---|---|
| White, n (%) | 41 (100) | 9 (100) | — |
| Married, n (%) | 36 (88) | 6 (67) | .14 |
| Tobacco use, n | 0 | 0 | — |
| Previous pregnancies, n a | 2.9 ± 1.8 | 3.4 ± 1.7 | .36 |
| ≥1 previous term delivery, n (%) | 17 (41.5) | 4 (44.4) | > .99 |
| Total number of previous preterm (20.0-36.9 wks’ gestation) deliveries, n a | 1.7 ± 0.8 | 2.1 ± 0.9 | .21 |
| History of preterm premature rupture of membranes, n (%) | 16 (39) | 2 (25) | .69 |
| Cervical cerclage placed in ≥1 pregnancies, n (%) | 11 (27) | 2 (22) | > .99 |
| Delivery gestational age of earliest previous preterm birth, wk a | 28.1 ± 4.0 | 31.5 ± 2.7 | .02 |
Responders delivered an average of +9.2 weeks later (range, +3.8 to +16.9 weeks) with 17P compared with +1.3 weeks (range, –1.9 to +2.9 weeks) for nonresponders ( P < .001). Two women delivered earlier with 17P: 1 woman delivered 11 days earlier, and the other woman delivered 13 days earlier. There were no indicated preterm deliveries in the study population. Pregnancy management was similar between responders and nonresponders. Nonresponders delivered earlier and were more likely to be preterm ( Table 2 ).
| Variable | 17-alpha hydroxyprogesterone caproate responders (n = 41) | 17-alpha hydroxyprogesterone caproate nonresponders (n = 9) | P value |
|---|---|---|---|
| Cervical cerclage placed, n (%) | 10 (24) | 2 (22) | > .99 |
| Vaginal cervical length assessed at least once during pregnancy, n (%) | 35 (92) | 8 (89) | > .99 |
| Delivery gestational age, wk a | 37.3 ± 1.8 | 32.9 ± 3.7 | < .001 |
| Delivered <37 weeks’ gestation, n (%) | 16 (39) | 8 (89) | .009 |
| Delivered <32 weeks’ gestation, n (%) | 0 | 3 (33) | .004 |
| Birthweight, g a | 3001 ± 544 | 2094 ± 755 | < .001 |
In our VAAST analysis, we allowed for recessive inheritance and locus heterogeneity and tested our genotypes using 1 million permutations. The genes with the greatest difference between responders and nonresponders are listed in Table 3 and represent the raw VAAST list of genes. The probability values that are displayed in Table 3 are unadjusted; none meet genome-wide significance ( P < 2.5 × 10 −6 ). Using the KEGG database as described earlier, we generated a list of 518 candidate genes in pathways that potentially are involved with SPTB and/or 17P metabolism ( Supplementary Table 1 ). This list of 518 genes was compared with the raw list of genes from our VAAST analysis. The NOS1 gene was the eighth highest scoring gene on the overall raw VAAST list and was the highest scoring variant from genes on the KEGG candidate gene list (VAAST, P < .00095).
| Rank | Gene | Gene name | P value a | VAAST | Gene function (if known) b |
|---|---|---|---|---|---|
| 1 | CUBN | Cubilin | 7.32e-5 | 23.3 | Receptor for intrinsic factor-vitamin B12 complexes |
| 2 | TMTC1 | Transmembrane and tetra-tricopeptide repeat containing 1 | 9.40e-5 | 17.1 | |
| 3 | TMEM158 | Transmembrane protein 158 | .00035 | 23.8 | Surface receptor proposed to function in neuronal survival pathway |
| 4 | ATMIN | ATM interactor | .00041 | 6.1 | |
| 5 | MYOG | Myogenin | .00056 | 13.0 | |
| 6 | NLRP10 | NLR family pyrin domain containing 10 | .00057 | 13.4 | Multifunctional negative regulator of inflammation and apoptosis |
| 7 | ENP1 | Essential nuclear protein 1 | .00094 | 13.8 | |
| 8 | NOS1 | Nitric oxide synthase 1 | .00094 | 11.2 | Synthesizes nitric oxide from L-arginine |
| 9 | PRMT6 | Protein arginine methyltransferase 6 | .0011 | 8.5 | Stimulates polymerase activity by enhancing DNA binding and processivity |
| 10 | ZFP28 | Zinc finger protein | .0023 | 10.7 | |
| 11 | CASZ1 | Castor zinc finger 1 | .00298 | 15.3 | Tumor suppression, blood pressure variation |
| 12 | VPS13C | Vacuolar protein sorting 13 homolog C | .00301 | 11.4 | |
| 13 | FCGR2A | Fc fragment of IgG, low affinity IIa receptor | .00305 | 14.6 | Found on phagocytic cells and involved with clearing of immune complexes |
| 14 | OR56A5 | Olfactory receptor family 56 subfamily A member 5 | .00318 | 13.0 | Smell perception |
| 15 | RPA4 | Replication protein A4 | .00337 | 5.5 | DNA double-strand break repair, inhibition of viral replication |
| 16 | ZNF853 | Zinc finger protein 853 | .00352 | 5.3 | |
| 17 | TSKS | Testis specific serine kinase substrate | .00373 | 8.0 | Tumorigenesis pathways and progression |
| 18 | MICAL2 | Microtubule associated monooxygenase, calponin and LIM domain containing 2 | .00384 | 10.2 | |
| 19 | CCDC50 | Coiled-coil domain containing 50 | .00401 | 11.1 | Hearing loss, effector of epidermal growth factor–mediated cell signaling |
| 20 | ANP32D | Acidic nuclear phosphoprotein 32 family member D | .00408 | 7.5 | Tumor suppressor |
| 21 | RALGAPA1 | Ral GTPase activating protein alpha subunit 1 | .00409 | 8.1 | |
| 22 | C16orf46 | Chromosome 16 open reading frame | .00409 | 5.9 | |
| 23 | ATP6V0A2 | ATPase lysosomal V0 subunit A2 | .00425 | 7.4 | Cutis laxa type II and wrinkly skin syndrome |
| 24 | SPRR1A | Small proline-rich protein 1A | .0457 | 4.9 | |
| 25 | PIH1D1 | PIH1domain containing 1 | .00464 | 2.1 |
a Unadjusted for multiple comparisons
b Gene functions per National Center for Biotechnology Information gene database ( ncbi.nlm.nih.gov ).
Next, the top 2.5% of raw genes (n = 457) that were generated from the initial VAAST analysis were selected for further analysis with the use of PANTHER ( Supplementary Table 2 ). These top 457 genes were compared with the online referent population that was provided by PANTHER. This analysis revealed several differentially expressed biologic processes ( Table 4 ) from our gene list compared with the referent population. The frequencies of genes that were classified by pathway or function in our gene list compared with expected proportions in a general population are also given in Table 4 . For example, based on gene distributions within the referent population, the expected frequency of genes in the nitric oxide pathway is 0.002*457 = approximately 1. However, from our list of 457 top genes that were identified by VAAST, the frequency of genes in the nitric oxide pathway was 0.02*457 = approximately 9 ( Table 4 ). We identified 8 genes in the nitric oxide synthase pathway with different allele frequencies between responders and nonresponders; our list included SPTA2, GA2L2, DMD, SYNE1, NOS1, MICA2, SMTL2, and DESP. All statistically significant pathways that are listed in Table 4 were over-represented in our gene list. There were no under- represented biologic processes from among our top gene list.
