Genetic effects on the timing of parturition and links to fetal birth weight

Phenotype definition

In this study, we included pregnancies with a singleton live birth and a spontaneous onset of delivery: medically initiated deliveries (either by induction or planned cesarean section) were excluded or part of controls for preterm delivery. Gestational duration in days was estimated using either the last menstrual period date or ultrasound. We excluded pregnancies lasting <140 days (20 completed weeks) or >310 days (44 completed weeks), as well as women with health complications prior to or during pregnancy and congenital fetal malformations. Spontaneous preterm delivery was defined as a spontaneous delivery <259 days (37 completed gestational weeks) or by using the ICD-10 O60 code, and controls as a delivery occurring between 273 and 294 days (39 and 42 gestational weeks). Post-term delivery was defined as a delivery occurring >294 days (42 completed weeks) or ICD-10 O48 code, and controls as a spontaneous delivery between 273 and 294 days (39 and 42 gestational weeks). Given the perfect genetic correlation between gestational duration and post-term delivery GWAS, and the small power of the latter, all downstream analyses are focused on gestational duration and preterm delivery.

Study cohorts and individual-level GWAS

This study consists of cohorts participating in the Early Growth Genetics (EGG) Consortium and the Norwegian Mother, Father and Child Cohort study (MoBa)³⁵, deCODE genetics¹⁰, Trøndelag Health Study (HUNT)³⁶, Danish Blood Donor Study (DBDS)³⁷, the Estonian Genome Center of the University of Tartu (EGCUT)³⁸ and summary statistics from FinnGen³⁹ and from a previous GWAS of gestational duration and preterm delivery performed using 23andMe data⁵. A total of 18 different cohorts (Supplementary Table 1) provided GWAS data under an additive model for meta-analysis for the maternal genome, resulting in 195,555 samples for gestational duration, 276,218 samples for preterm delivery (n cases = 18,797) and 131,279 samples for post-term delivery (n cases = 15,972) of recent European ancestries (indicated by principal component analysis). For binary outcomes (preterm and post-term deliveries), only cohorts with an effective sample size >100 were included. Detailed description of the cohorts included can be found in the Supplementary Note. All study participants provided a signed informed consent, and all research studies were approved by the relevant institutional ethics review boards (Supplementary Note).

Each individual cohort applied specific QC procedures, data imputation and analysis independently following the consortium recommendations. Unless more stringent, samples were excluded if genotype call rate <95%, autosomal mean heterozygosity >3 standard deviations from the cohort mean, sex mismatch or major recent ancestry was other than European (HapMap central European). Genetic variants were excluded if genotype call rate <98%, Hardy-Weinberg equilibrium P value < 1 × 10⁻⁶ or minor allele frequency <1%. Reference panels for imputation were either 1000 Genomes Project⁴⁰, Haplotype Reference Consortium⁴¹, 10KUK or a combination of one of the mentioned reference panels and own whole-genome sequencing data (deCODE, HUNT, DBDS and FinnGen). Each individual cohort performed a GWAS using an additive linear regression model adjusted for, at least, genetic principal components or relationship matrix on autosomal chromosomes and chromosome X. Summary statistics for each individual cohort were stored centrally and underwent QC procedures before meta-analysis (Supplementary Note).

Meta-analysis of GWASs

After QC, individual-cohort GWAS summary statistics were pooled using fixed-effects inverse-variance weighted meta-analysis with METAL⁴² without genomic control correction. We also performed an analysis of heterogeneity of effects (Supplementary Table 2; I² statistic). After meta-analysis, we removed genetic variants reported in less than half the number of available samples for each phenotype, resulting in 9-10 million genetic variants. For example, the variant observed in the largest number of samples for gestational duration was available in 195,555 individuals; only variants reported in at least 97,778 were kept. Genomic inflation factors were low for all three phenotypes (Supplementary Table 16; gestational duration λ = 1.14, preterm delivery λ = 1.08 and post-term delivery λ = 1.05). LD-score regression intercepts were substantially lower than genomic inflation factors, suggesting that the inflation in test statistics was mostly due to polygenicity (Supplementary Table 16). Test statistics were not further adjusted for genomic control for any of the phenotypes. If not otherwise stated, all analyses presented in this study are two-sided tests.

Initially, we naively defined independent loci based on physical distance, where SNPs within 250 kb from the index SNP were considered to be at the same locus. Novel loci were defined as loci not overlapping previously reported gestational duration loci in the largest GWAS performed to date⁵. Finally, we used conditional analysis to resolve independent loci (see below).

Conditional analysis

We looked for conditionally independent associations within each locus using approximate conditional and joint (COJO) analysis⁴³ implemented in Genome-wide Complex Trait Analysis (GCTA) software⁴⁴. We ran a stepwise model selection (-cojo-slct) to identify conditionally independent genetic variants at P < 5 × 10⁻⁸ for each of the genome-wide significant loci (using a radius of 1.5 Mb from the index SNP). Overlapping loci were merged into a single locus (only two loci overlapped, at 3q23). LD between genetic variants was estimated from 19,092 maternal samples from the Norwegian Mother, Father and Child Cohort, after excluding variants with imputation INFO score <0.4. We converted the reference panel from BGEN files to hard-called PLINK binary format (.bed). As per default in COJO, genetic variants >10 Mb apart were assumed to be in complete linkage equilibrium.

Gene prioritization

To prioritize genes at the gestational duration loci identified, we set the baseline as the nearest protein-coding gene to the index SNP at each independent locus. Although naive, this approach has been consistently shown to outperform other single metrics for locus-to-gene mapping^45,46. Next, we performed colocalization analysis for cis-eQTLs in 1,367 human induced pluripotent stem cell lines from the i2QTL resource (±250 kb from gene start and stop position)¹⁵, endometrium (± 250 kb from gene start and stop position)¹⁶ and uterus, vagina and ovary from GTEx (±1 Mb around transcription start site)¹⁷. None of the variants we identified were in LD (r² > 0.6) with missense variants. To complement the prioritization of genes, we queried each of the index SNPs for blood protein QTLs¹⁸ (both in cis and trans). For all index SNPs that were protein QTLs (P < 5 × 10⁻⁶), we performed colocalization analyses (±1.5 Mb around the index SNP). We excluded the HLA region due to its large pleiotropic effects.

Colocalization

We utilized genetic colocalization to identify pleiotropic effects between gestational duration and expression and protein quantitative trait locus (see Gene prioritization) and with other female and reproductive traits. To this end, we applied COLOC¹⁴, which evaluates, in a Bayesian statistical framework, whether a single locus from two different phenotypes best fits a model where the associations are due to a single shared variant or distinct variants in close LD (Supplementary Note).

Prior probabilities for each for the non-null hypotheses were set as suggested by Wallace (prior probabilities that a random SNP in the loci is associated with phenotype A, phenotype B, or both phenotypes, 1 × 10⁻⁴, 1 × 10⁻⁴, and 5 × 10⁻⁶, respectively), which are considered more conservative than the ones set by default⁴⁷. Strong evidence of colocalization was defined as a posterior probability of colocalization >0.9.

Enrichment analysis

We tested for enrichment based on top loci and genome-wide using partitioned LD-score regression. To test for overrepresentation in tissue-specific RNA expression (Human Protein Atlas, RNA consensus tissue gene data)⁴⁸, a Wilcoxon rank-sum test was performed on normalized RNA for genes within our set (above-mentioned) and all other genes. Significance for this test was set at Bonferroni correction for the number of tissues (P < 0.05/61), and suggestive evidence at P < 0.1/61. At the genome-wide level, we performed partitioned heritability using LD-score regression to test for enrichment in 97 different annotations^49,50, tissue-specific RNA expression using 205 different tissues/cell types⁵¹, using precomputed partitioned LD-scores for subjects of recent European ancestry (baseline-LD model v2.2) and for enrichment in regions harboring genes differentially expressed during labor in single cells from myometrium²².

Genetic correlations

We estimated genetic correlations by performing LD-score regression⁵² locally using precomputed LD-scores from 1000 Genomes Project samples of recent European ancestry. The MHC region (chr6:28477797-33448354) was removed prior to running LD-score regression.

Resolving effect origin

To classify the identified index SNPs for gestational duration as having maternal, fetal, or maternal and fetal origin, we performed an association analysis using the parental transmitted and nontransmitted alleles on gestational duration. We used phased genotype data (that is estimated haplotypes) in parent-offsprings or mother-child duos to infer the parent-of-origin of the genotyped/imputed alleles as previously described³⁰. Once the transmitted allele was identified, the nontransmitted maternal allele was extracted. Briefly, parental origin of each allele was inferred using genotypes of relatives, reference cohort data, or distributions of genotypes within the cohort and LD measurements. Different methods were used for phasing in each of the cohorts providing data for this analysis^{10,53,54,55,56} (Supplementary Table 17). Details of the phasing strategy used by each cohort are described in Supplementary Note.

For each index SNP, we fit the following linear regression model:

where MnT and MT refer to the maternal nontransmitted and transmitted alleles respectively, and PT refers to the paternal transmitted alleles. The latter is interpreted as a fetal-only genetic effect, whereas the effect of the maternal nontransmitted allele is a maternal-only genetic effect. We first estimated the effects of the index SNPs in each birth cohort separately; effect sizes were then combined through fixed-effect meta-analysis, totaling a sample size of 136,833 (Supplementary Note and Supplementary Table 17). To classify the identified genetic variants into classes with similar patterns of effect, we used model-based clustering¹⁰. Variants were clustered based on estimated effects of the transmitted and nontransmitted parental alleles into five clusters. Two clusters assume fetal effect only, one with effect independent of parent of origin, and one where the effect is limited to the maternally transmitted allele; a cluster with maternal effect only; and two clusters with both maternal and fetal effects, either in opposite or same direction.

Locus pleiotropy at 3q21

After identifying locus pleiotropy between the maternal effect on gestational duration and the fetal-only effect on birth weight at the ADCY5 gene region, we set out to investigate differences between the two top SNPs in their colocalization with other traits. Phenome-wide colocalization for the two regions was performed using summary statistics from FinnGen (data freeze 5) and Pan UK Biobank data (https://pan.ukbb.broadinstitute.org, in participants of recent European ancestry; Supplementary Note).

Female reproductive traits

We obtained summary statistics for several female reproductive traits from different sources (minimum sample size 10,000). We included summary statistics from the following traits: miscarriage²⁶, gestational duration (fetal genome)⁶, age at first birth, age at menarche (Neale lab, http://www.nealelab.is), age at menopause⁵⁷, number of live births (Neale lab, http://www.nealelab.is), testosterone⁵⁸, CBAT⁵⁸, SHBG⁵⁸, estradiol (women, Neale lab, http://www.nealelab.is), pelvic organ prolapse (FinnGen), polycystic ovary syndrome (⁵⁹ and FinnGen), endometriosis (Neale lab, http://www.nealelab.is), leiomyoma uterus (FinnGen) and pre-eclampsia⁶⁰. For polycystic ovary syndrome, we meta-analyzed summary statistics from the largest published GWAS⁵⁹ and FinnGen. We estimated genetic correlations between gestational duration and preterm delivery and these traits, and latent causal variable analysis between sex hormones (testosterone, CBAT and SHBG) and gestational duration and preterm delivery. We further explored causality using two-sample Mendelian randomization and inspected whether the effects originated in the maternal or the fetal genome (see below, ‘Mendelian randomization’). Finally, when one trait is causally upstream of the other, it is expected that the two traits would share a causal variant at some of the trait-associated loci. To test for this at the genome-wide scale, we performed colocalization analysis between sex hormones and gestational duration and preterm delivery using approximately LD-independent regions⁶¹.

Gestational duration and preterm delivery polygenic scores

To obtain an independent sample for training and validation of a polygenic score, the meta-analyses for gestational duration were rerun, excluding the MoBa cohort. These new meta-analysis results were used as the base data sets to calculate the polygenic scores. After applying the same exclusion criteria as used for the study samples in the meta-analysis, and removing duplicated samples and those with a kinship of greater than 0.125, the MoBa cohort was randomly split, using 80% (n = 15,768) as the training cohort and the remaining 20% (n = 3,942) as the validation cohort. LDpred2 was used for the calculation of the polygenic scores²⁵. A description of polygenic score training can be found in Supplementary Note.

Polygenic score validation

We constructed polygenic scores converted to z-scores to enable comparison of the gestational duration and the preterm delivery polygenic scores. To test the performance of the polygenic score, a linear regression was conducted for gestational duration by the polygenic score. A second model was used that adjusted for five principal components and genotyped batch. R² was calculated for the models to quantify variance explained.

The utility of the polygenic score for the prediction of preterm delivery was also assessed. Gestational duration was dichotomized into preterm delivery (<37 weeks) or full term (≥39 weeks and <41 weeks). Two models were analyzed, one assessing just the polygenic score and a second adjusting for five principal components and genotype batch. Receiver operating characteristic, area under the curve were calculated for each model and used as assessment of diagnostic accuracy.

Mendelian randomization

We performed Mendelian randomization to study the effects of gestational duration (maternal) on birth weight (maternal) and the effects of fetal growth (fetal effect on birth weight) and sex hormones on gestational duration.

To study the effect of gestational duration on birth weight, we employed two-sample Mendelian randomization. The 24 index SNPs (22 autosomal SNPs) from the present gestational duration meta-analysis and the effect sizes from the parental transmitted and nontransmitted alleles analysis were used to instrument gestational duration. Birth weight was instrumented using summary statistics from a previous GWAS of offspring’s birth weight with minimal adjustment by gestational duration (<15% of samples)⁹.

We assessed the effect of sex hormones (testosterone, SHBG and CBAT) on gestational duration using two-sample Mendelian randomization and instrumenting the hormones using a polygenic score for the parental transmitted and nontransmitted alleles. For each sex hormone, we obtained a list of independent SNPs genome-wide associated with these traits (Supplementary Table 9) by performing GWAS clumping (r² > 0.001) using the following PLINK command:

plink–bfile <1000 Genomes > –clump {GWAS summary statistics}–clump-r2 0.001–clump-kb 1000–clump-p1 5e-8–clump-p2 1e-5.

We also used a set of SNPs associated with testosterone, but with no aggregated effects on SHBG, as clustered in²⁸. Such variants were used as instrumental variables in the two-sample Mendelian randomization analysis and to construct the polygenic score for the parental transmitted and nontransmitted alleles. The current meta-analysis results were employed as outcome for the two-sample Mendelian randomization analysis (inverse-variance weighted and MR-Egger). We subsequently constructed the polygenic score for the maternal transmitted and nontransmitted alleles and the paternal transmitted alleles in 46,105 parent-offsprings from Iceland and Norway. We estimated the effects of the maternal nontransmitted (MnT_PGS) and transmitted (MT_PGS) and paternal transmitted (PT_PGS) alleles polygenic score using the following linear model:

Again, effects from each of the three data sets (Iceland, MoBa and HUNT) were combined using fixed-effect inverse-variance weighted meta-analysis.

To understand the impact of fetal growth on gestational duration, we used individual genetic data from 35,280 (ultrasound-gestational duration) and 48,741 (last menstrual period-gestational duration) parent-offsprings from Iceland, the MoBa cohort and HUNT. To instrument fetal growth, we used 68 SNPs with fetal-only effect on birth weight as classified in Warrington et al. using Structural Equation Modeling⁹. Based on these 68 SNPs, we constructed a fetal growth polygenic score for the parental transmitted and nontransmitted alleles and regressed these on gestational duration (estimated by ultrasound or last menstrual period, separately). We estimated the effects of the maternal nontransmitted (MnT_PGS) and transmitted (MT_PGS) and paternal transmitted (PT_PGS) alleles polygenic scores as above.

Effect estimates from each of the three data sets (Iceland, MoBa and HUNT) were pooled using fixed-effects inverse-variance weighted meta-analysis.

Multitrait conditional analysis

GCTA was used to perform bi-directional multitrait COJO (mtCOJO)²⁹ analysis using summary statistics. The gestational duration GWAS was conditioned on the birth weight GWAS and vice versa (Supplementary Note), using birth weight summary statistics from the largest GWAS meta-analysis of birth weight⁹. We did not condition on the fetal effects on gestational duration due to a lack of power in the fetal GWAS⁶.

Maternal–fetal pleiotropy on gestational duration and birth weight

We further investigated what are the fetal effects on birth weight for the maternal gestational duration-increasing alleles, and the maternal effects on gestational duration for the fetal birth weight-increasing alleles. To study this, we borrowed the inverse-variance weighted analysis from Mendelian randomization, but using the effects of two distinct genomes, the maternal and fetal. We caution that this should not be interpreted under a causal framework.

To understand what the maternal gestational duration-raising alleles do to birth weight when present in the fetus, we used the effect sizes and standard errors of the parental transmitted and nontransmitted alleles for the 22 autosomal index SNPs on gestational duration and assessed its effects on the same SNPs with a fetal-only effect on birth weight. To understand what the fetal birth weight-raising alleles do to gestational duration when present in the mother, we used the effect sizes and standard errors of 68 autosomal SNPs associated with fetal effects on birth weight and the effect sizes and standard errors from the current maternal GWAS of gestational duration.

Evolutionary analysis

To examine the evolutionary history of the regions identified in the GWAS meta-analysis, we ran the significant variants through the MOSAIc pipeline²³. This pipeline is designed to detect enrichment in evolutionary signals using a variety of sequence-based metrics of selection (Supplementary Note).

Variant annotation

Variants were annotated using Ensembl’s Variant Effect Predictor (hg19) command line tool⁶². Physical coordinates of protein-coding genes were obtained from the UCSC Table Browser⁶³, and were matched to the index SNPs using bedtools v2.29.2 (ref. ⁶⁴).

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.