Genetic Polymorphisms Associated with Adverse Pregnancy Outcomes in Nulliparas

Participants

The participants of the analysis were enrolled in the nuMoM2b cohort (https://www.nichd.nih.gov/research/supported/nuMoM2b), a longitudinal, multiethnic cohort study of nulliparous individuals. All participating centers obtained approval by the local Institutional Review Boards (IRBs) of their corresponding recruitment institutions, with further details on data collection having been previously described^11,12. The study enrolled 10,038 nulliparous people from the first trimester of their pregnancy to participate in three study visits during pregnancy. Collection of health status and biomarkers were conducted at regular intervals, and documentation of pregnancy outcomes was performed by chart abstraction using a priori definitions. The details of this process were described elsewhere⁷.

Phenotype definitions

Gestational length: We opted to use a quantitative phenotype, gestational length, instead of a binary preterm/full term outcome to gain additional information and statistical power by using a more granular phenotype. Gestational length was determined from an estimated due date established by a first-trimester ultrasound crown-rump length measurement and was recorded in weeks¹². Preterm birth was defined as any live birth that occurred before 37 weeks gestational age. Cases of stillbirth, fetal demise, elective termination, and indicated termination were all excluded from this phenotype group (Table 1).

Preeclampsia: Cases include diagnosis of preeclampsia (with and without severe features), eclampsia, and chronic hypertension with superimposed preeclampsia. A detailed description of nuMoM2b study definitions of hypertensive disorders of pregnancy was published in the supplement to the paper by Facco et al¹³. All other individuals were treated as controls.

Gestational diabetes (GDM): GDM was diagnosed through clinical evaluation, from fasting blood sugar, sequential 1-hour glucose challenge test followed by a 3-hour glucose tolerance test (GTT), or a single step 2-hour 75-gram GTT¹³. We excluded individuals diagnosed with pregestational diabetes from GDM analyses, and all other individuals were treated as controls. Pregnancy loss: All subjects who had a pregnancy loss, regardless of gestational age, were considered as cases (Table 1). Individuals who underwent termination of pregnancy were excluded from the pregnancy loss analysis. Subjects who had a live birth were treated as controls.

Genotyping

The genotyped cohort comprised 9,757 participants who had adequate samples and agreed to be genotyped. DNA extractions from whole blood were carried out on a QIAsymphony instrument (from Qiagen; extraction kit DSP DNA Midi Kit #937355, protocol Blood_1000_V7_DSP) at the Center for Genomics and Bioinformatics (Indiana University, Bloomington). Genotyping was completed at the Van Andel Institute (Grand Rapids, MI, USA) using the Infinium Multi-Ethnic Global D2 BeadChip (Illumina, Miami, USA). We imposed standard filters for quality control of loci at this stage (cluster separation < 0.3, AA R Mean < 0.2, AB R Mean < 0.2, BB R Mean <0.2, 10% GC < 0.3) using GenomeStudio v2.4 (Illumina). Genotype calls (in GCT format) for the 1,748,280 loci that passed initial quality control were made with Beeline autoconvert (Ilumina).

Inclusion criteria, QC pipeline, Ancestry determination, and imputation

We performed preprocessing to determine any sex inconsistencies, autosome missingness (>5%), and contamination (more than 25 related subjects within the same dataset). Given the highly heterogeneous data set, we used the KING-Robust algorithm, a pairwise kinship estimator for GWAS that is robust to the presence of unknown population substructure¹⁴. We inferred kinship estimates between all pairs of subjects, randomly removing one subject from the pairs of subjects with first- or second-degree relatedness such that we minimized the number of subjects removed. The remaining samples were filtered to remove SNPs with a minor allele frequency (MAF) less than 0.01, genotyping rate <95%, and a Hardy-Weinberg Equilibrium (HWE) P < 5×10⁻⁶.

We determined the ancestry of each subject using SNPweights v.2.1 and leveraging approximately 40,000 ancestry informative markers available from the 1000 Genomes Consortium¹⁵. Applying a probability cut-off of >50%, the samples were clustered into five ancestry groups concordant with the 1000 Genomes Consortium¹⁵ super-populations: African (AFR, n=1425), American (AMR; n=846), East Asian (EAS; n=323), European (EUR, n=6082), and South Asian (SAS; n=112). We observed a large fraction of highly admixed individuals, thus we also established a sixth group (ADM, n=891) of subjects who have none of the five ancestry group percentages reaching >50% (Supplemental Figure S1). The heterogeneous genetic ancestry of this cohort would result in high levels of genomic inflation in a single GWAS analysis. To reduce this confounding effect, we clustered the cohort by genetic ancestry for downstream analyses (Supplemental Figure S1, Supplemental Figure S2). Thus, from this point subjects were analyzed in separate groups according to their majority genetic ancestry when doing the filtering preprocessing and imputation.

Not all the ancestry groups contained enough subjects to power a GWAS study. Thus, we only proceeded with imputing the European, African, and American ancestry subjects. These subjects were phased with EAGLE and imputed by Minimac3¹⁶ using the TOPMED Imputation Server^17,18. We used version R2 of the TOPMED panel, currently the largest panel of sequenced human genomes, and containing all the ancestries we observed in our cohort. Following imputation, we excluded SNPs with an r² quality score <0.7, MAF <0.01, genotyping rate <95%, and a Hardy-Weinberg Equilibrium (HWE) P < 5×10⁻⁶ within each imputed group.

Population structure was determined by PLINK using pruned SNPs from the data (Linkage Disequilibrium pruning r² < 0.1). The top ten principal components were computed for each of the three sub-cohorts (EUR, AFR, AMR). As maternal age has a non-linear effect on preterm birth, we transformed each subject’s age into the distance from the median percentile in all subjects.

Genome-wide associations

Association testing was carried out using regression models implemented in PLINK v1.9¹⁹. The suggestive association threshold was P < 1×10⁻⁵ and the threshold for genome-wide significance was P < 5×10⁻⁸. The model was adjusted for each subject’s rank-transformed age and the first ten principal components from our population structure analysis. We used all ten computed principal components due to the high level of population substructure and genetic admixture present in the nuMoM2b data set. Any subjects with missing age or phenotype information were excluded from the analysis. Results from each sub-cohort were fixed-effect meta-analyzed using GWAMA²⁰ and SNPs that did not appear in at least two of the three individual GWAS were removed. Result plots were displayed using R libraries. The SNPs of interest were further annotated with Variant Effect Predictor²¹ to identify their corresponding gene function. The model for gestational length was adjusted for rank-transformed age.

Using the results of the genome-wide associations for each APO, we inferred SNP-based heritability of these traits using LDscore²². We searched for phenotypes that had significant genetic correlations with the APOs using the LD-Hub database²³. Finally, we scanned for putative regulatory effects by carrying out a transcriptome-wide association analysis (TWAS) as implemented in the Fusion software, and using the available expression reference weights for whole blood and adipose tissue²⁴, as well as liver, pancreas, vagina, and uterus, and whole blood from the GTEx v7 multi-tissue RNA-seq data set. For TWAS, P-values were Bonferroni-corrected by the number of genes in each panel. Conditional expression analyses were performed using R scripts from the Fusion software.

Pregnancy Loss

This trait showed the highest SNP-based heritability (0.25) of the four APOs studied. We observed two novel associations with pregnancy loss (Table 2): rs62021480 and rs142795512. Located on chromosome 15, rs62021480 encodes a synonymous variant on RMGA (Repulsive Guidance Molecule BMP Co-Receptor A), a gene that encodes a glycoprotein that guides developing axons and may act as a tumor suppressor. Rs142795512 is an intronic variant in the gene TRMU (TRNA 5-Methylaminomethyl-2-Thiouridylate Methyltransferase), which encodes a mitochondrial tRNA-modifying enzyme. Variants in TRMU have been linked to several disease phenotypes, including infantile liver failure²⁵.

Through TWAS, we found a correlation between gene expression levels of TTC38 (Tetratricopeptide Repeat Domain 38) in uterine tissue and pregnancy loss (Z=4.83, P=1.36×10⁻⁶; Supplemental Figure S10).

Gestational Length

The genomic inflations of each sub-cohort GWAS ranged between λ = 0.99 and λ = 1.02, while the GWA meta-analysis across sub-cohorts did not show genomic inflation (Supplemental Figure S3). SNP-based heritability for this trait was 23%. A single SNP, rs144380855, reached genome-wide significance in the meta-analysis of gestational length. rs144380855 appears to be associated with a reduced gestational length (beta = -2.07, P = 5.3×10⁻¹⁰), and is an intergenic SNP (Table 2). We found evidence of association (P < 0.05) in 8 out of 34 SNPs previously reported to be associated with gestational length (Table 3; Supplemental Table S1). These 8 markers were clustered in five genes: WNT4 (Wnt Family Member 4), EBF1 (EBF transcription factor 1), PER3 (Period Circadian Regulator 3), IL10 (Interleukin 10), and ADCY5 (Adenalyte Cyclase 5).

Gestational Diabetes mellitus

Three significant loci were identified in the GWAS for GDM (Table 2). The top locus included intronic variants rs61167087 and rs72956265, in the gene ZBTB20, a Zinc Finger and BTB Domain Containing 20. The second significant association is also an intronic SNP, rs117689036 in MDGA2, a MAM Domain-Containing Glycosylphosphatidylinositol Anchor 2. Finally, rs10890563 is a 3’UTR variant on the gene GUCY1A2, Guanylate Cyclase 1 Soluble Subunit Alpha 2. Several loci show suggestive associations in the meta-analysis, which we report in the supplemental text (Figure 4; Supplemental Table S2). We did not find genomic inflation in any of the ancestry groups or in the meta-analysis (λ =1.004, Supplemental Figure S5), and we estimated SNP-based heritability to be 0.17.

We also found a positive association consistent with an effect on gestational diabetes for one previously reported variant (rs4506565, P = 0.001; Table 3). However, we failed to find any confirmatory signal in 17 other variants (Supplemental Table S1).

Preeclampsia

We found no SNPs associated with preeclampsia (Supplemental Figure S6). Again, genomic inflation was low (λ =1.005, Supplemental Figure S7), and our SNP-based estimate of heritability (0.02) suggests a weak genetic basis for this APO. We confirmed 1 of 22 previously linked variants: rs17367504 (near gene MTHFR) showed a consistent effect across our three studied sub-populations (Table 3; Supplemental Table S1).

Principal Findings

We identified a novel locus associated with gestational length, along with three loci associated with GDM and two additional loci associated with the previously understudied phenotype of pregnancy loss. Our study suggests that variants in ZBTB20, a gene that has been previously found to regulate beta cell function and glucose homeostasis in mice²⁶, may be involved in the development of GDM. Furthermore, TWAS analysis implicates the gene TTC38 in pregnancy loss across four tissue types (liver, ovary, pancreas, and uterus). Consistent with this, one of the top SNPs in the pregnancy loss GWAS (rs114058457, Supplementary Table S2), was previously established by Jansen et al. as an eQTL for TTC38 in blood in a cohort of 4,896 subjects of European ancestry (FDR-corrected P < 1.34×10⁻⁵).²⁷

Clinical Implications

The SNPs that reached genome-wide significance in our study highlight the role of previously unknown genes in our APOs of interest. These genes may further implicate specific biological interactions and pathways that are disrupted in each APO. While not reaching genome-wide significance, we were able to document the directionality of association for several SNPs previously reported for other APOs. These SNPs, now having been discerned by different investigations, are valuable by identifying genes relevant to pathophysiology–which can lead to better prevention and treatment.

Clinical utility of genetic screening is most effective when there is a preventive measure that can be employed. If strong genetic markers for preeclampsia, for instance, had been discovered in pathway genes for mechanisms that could be influenced by aspirin therapy, these would be prime candidates for screening nulliparas who might not have other risk factors that would lead to recommended prophylaxis ²⁸. SNPs implicated in shorter gestational length could lead to further screening with cervical length ultrasound or perhaps the administration of preventive progestin therapy. Incorporating genetic screening for APOs into clinical care will necessitate more data and trials, but novel and confirmatory findings such as those presented here are initial steps towards potential clinical utility.

Research Implications

Given the low heritability computed across all four APOs, it appears that environmental or other factors may play an important role in these phenotypes. One possibility is that individual genetic loci each confer small effects to the phenotype as a whole; indeed, our findings are consistent with this hypothesis. The development of statistical methods that leverage such assumptions, such as polygenic risk scoring, may enable better predictive capacities of models. Evolutionarily, given that variants with a strong effect on infant mortality may face strong negative selection, it may be more likely for rare variants of smaller effect to play a larger role in the genetics of preterm birth²⁹. Identifying rare causative variants would require sequencing efforts including whole-genome and whole-exome sequencing, as well as larger study populations of diverse genetic ancestry. Additionally, our study has unearthed potential contributors to preterm birth at multiple levels of granularity—SNP, gene, and tissue. Targeted follow-up functional studies of these layers of genetic factors in vivo or in vitro may better elucidate the underlying biological mechanisms influencing preterm birth phenotypes and ultimately can shape targeted treatments.

Strengths and Limitations

While ours is one of the largest GWAS on traits associated with pregnancy (e.g., it is, to our knowledge, the largest of three GWAS on GDM^30,31), it highlights the limitations of prospective studies in a population of highly diverse genetic ancestry. Namely, this is an imbalanced cohort both in terms of cases-controls and genetic ancestries. Our results would benefit from further validation cohorts and deeper sampling of cases, particularly of individuals of non-European genetic ancestry. Larger sample sizes and future meta-analysis would allow us to explore the diversity of types and etiologies of preterm birth and other APOs. As we noted previously, the environment may play a large role in mediating our APOs, and therefore future directions include conducting gene-by-environment interaction studies to further elucidate the interplay between the two.

Conclusions

Our study identified novel genetic associations for three pregnancy phenotypes, reporting new variants correlated with gestational length, gestational diabetes, and pregnancy loss. This study confirms the presence of genetic signals in previously identified pathways for preterm birth, GDM, and preeclampsia. Validation of these new associated SNPs in other cohorts will be important in determining their potential clinical usefulness.