Loss of Imprinting in Human Placentas Is Widespread, Coordinated, and Predicts Birth Phenotypes.

Genomic imprinting leads to mono-allelic expression of genes based on parent of origin. Therian mammals and angiosperms evolved this mechanism in nutritive tissues, the placenta, and endosperm, where maternal and paternal genomes are in conflict with respect to resource allocation. We used RNA-seq to analyze allelic bias in the expression of 91 known imprinted genes in term human placentas from a prospective cohort study in Mali. A large fraction of the imprinted exons (39%) deviated from mono-allelic expression. Loss of imprinting (LOI) occurred in genes with either maternal or paternal expression bias, albeit more frequently in the former. We characterized LOI using binomial generalized linear mixed models. Variation in LOI was predominantly at the gene as opposed to the exon level, consistent with a single promoter driving the expression of most exons in a gene. Some genes were less prone to LOI than others, particularly lncRNA genes were rarely expressed from the repressed allele. Further, some individuals had more LOI than others and, within a person, the expression bias of maternally and paternally imprinted genes was correlated. We hypothesize that trans-acting maternal effect genes mediate correlated LOI and provide the mother with an additional lever to control fetal growth by extending her influence to LOI of the paternally imprinted genes. Limited evidence exists to support associations between LOI and offspring phenotypes. We show that birth length and placental weight were associated with allelic bias, making this the first comprehensive report of an association between LOI and a birth phenotype.

Introduction

Transcription of imprinted genes is repressed on either the paternal or the maternal allele, resulting in allele-specific expression (ASE). Imprinting evolved in therian mammals and angiosperms, both of which have embryo nourishing tissues, the placenta, and the endosperm, respectively. In humans, a small fraction of genes (∼100) are imprinted; many are imprinted in tissues throughout the body while others are exclusively imprinted in the placenta (Peters 2014). Conservation of imprinting among mammalian species is incomplete as evidenced by genes that are only imprinted in some taxonomic groups, such as primates (Noguer-Dance et al. 2010).

The kinship hypothesis explains the evolutionary forces that led to imprinting and its maintenance (Moore and Haig 1991). In essence, it attributes imprinting to a conflict of interest between the maternal and paternal genomes over resource allocation to the fetus. This hypothesis successfully predicted that growth promoting genes are usually repressed on the maternal allele while growth inhibiting genes are repressed on the paternal allele (Babak et al. 2015); it also explained why the phenotypes of uniparental conceptuses (McGrath and Solter 1984; Surani et al. 1984) are excessively small or large, depending on the parent. A growing body of experimental validation for the kinship hypothesis is accumulating from the growth phenotypes of knock-out mice (Fowden et al. 2011) and RNA-seq expression analysis (Babak et al. 2015). When homologous genes are compared, those that are imprinted in mammals are expressed at higher levels from one allele than from two alleles in species lacking imprinting (chicken and platypus), a finding that has provided further evidence for positive selection and conflict between maternal and paternal genomes (Babak et al. 2015).

Historically, the imprinting status of a gene was categorized as either paternal or maternal based on the analysis of a few informative SNPs with semiquantitative methods. Such early assays already revealed heterogeneous imprinting in multiple tissues as well as variability between individuals (Dao et al. 1998; Sakatani et al. 2001; McMinn, Wei, Sadovsky et al. 2006). Recently, RNA-seq has enabled the genome-wide assessment of allelic bias in human tissues (Metsalu et al. 2014; Babak et al. 2015; Baran et al. 2015; Hamada et al. 2016; Gulyas-Kovacs et al. 2018; Zink et al. 2018; Jadhav et al. 2019) and has shown that repression is often incomplete with some expression from the silenced allele. In complex tissues, specious loss of imprinting (LOI) might arise if only a fraction of cell types express the gene from one allele while other cell types in the tissue have biallelic expression. However, LOI is also detected in single cell RNA-seq studies demonstrating that LOI can be due to derepression (Santoni et al. 2017; Vertesy et al. 2018).

Two studies aimed at the discovery of imprinted genes by RNA-seq analysis in various tissues have produced largely overlapping gene lists in humans (Babak et al. 2015; Baran et al. 2015) and the same is true for mice (Babak et al. 2015; Bonthuis et al. 2015; Perez et al. 2015) suggesting that the majority of imprinted genes have been identified in most tissues including the placenta (Metsalu et al. 2014; Hamada et al. 2016). The mechanistic understanding is less complete even though it is well established that differentially methylated regions (DMR) are essential for repression at many loci (Bartolomei and Ferguson-Smith 2011). However, additional imprints, like H3K27me3, have also been shown to be functional (Inoue et al. 2017). Since the expression of imprinted genes is influenced by multiple mechanisms, analysis of RNA provides a more inclusive readout than DNA methylation to characterize the functional variation of imprinted genes.

Here, we report the first systematic analysis of the variation in ASE of imprinted genes in the placenta for a human population. We generated RNA-seq data from placentas collected at the time of birth from participants in a multigenerational cohort in Mali, West Africa. We computed allele-specific read counts for all eSNPs that mapped to known imprinted genes. The resulting data from 52 placentas were analyzed using binomial generalized linear mixed models (GLMMs) using the imprinting call at the exon level as the dependent variable.

We use the term LOI to refer to departures from mono-allelic RNA expression as determined by RNA-seq and pyrosequencing rather than to measures of the loss of molecular imprints, such as DNA methylation. LOI occurred in both maternally and paternally silenced genes, but was more pronounced in paternally silenced genes. LOI was not uniform across genes, with some genes prone to LOI (RHOBTB3, TFPI2) and others rarely showing LOI (KCNQ1OT1, MEG3). Imprinted lncRNA genes showed less expression from the silenced allele than did protein coding genes. Unexpectedly, we found a significant correlation in LOI between paternally expressed genes (PEGs) and maternally expressed genes (MEGs) in the same placenta. Some individuals had more LOI than others, consistent with a role for LOI in shaping phenotypic outcomes. Birth length and placental weight were associated with allelic bias, making this the first report to analyze interperson variation in LOI across imprinted genes in humans and to find an association with birth phenotype. We conclude that LOI is a nonstochastic process with significant gene-level variation as well as person-level variation that has phenotypic consequences.

New Approaches

Sample Collection and Sequencing

We collected placental specimens from participants in a multigenerational, population-based cohort study of the Dogon of Mali, West Africa (see Materials and Methods). When the young women (F1) who had been followed throughout childhood and adolescence gave birth, placental and fetal cord blood samples were collected and the infant was enrolled as a member of the F2 generation. Here, we report data for 52 placentas collected from 52 mothers as well as birth parameters for the F2 infants. Descriptive statistics are summarized in supplementary table S1, Supplementary Material online.

The placentas were collected by trained midwives after natural delivery in a local hospital. Compared with placental cohorts from high income countries (Kappil et al. 2015; Moore et al. 2015; Prats-Puig et al. 2017), the mothers in this study were sub-Saharan Africans who shared the same ethnicity, were younger (mean age 18.2 years ± 1.4), predominantly primiparous (85%), and the mean weight of their newborns (2.6 kg ± 0.4) and their placentas was lower (475 g ± 87). We sampled two cotyledons from each placenta and analyzed the samples separately.

To accommodate partially degraded RNA, we employed stranded sequence library methods with random priming and we sequenced to high depth (∼300×106 reads, see supplementary table S2, Supplementary Material online and Materials and Methods).

Identification of hetSNPs and Assignment of Imprinting Codes

As most imprinted genes in humans have been identified, we were able to design a targeting region for genotyping that included 91 genes known to be imprinted and to assign the parental origin of the imprint. In the targeted region, we located all heterozygous eSNPs from DNA isolated from umbilical cord, which is a fetal tissue free of maternal DNA. We coded genes based on the allele expressed preferentially as Paternal (PEGs) (N = 40), Maternal (MEGs) (N = 22), or Complex (CEGs) (N = 29); CEGs had gene isoforms with dissimilar parental expression bias or conflicting reports in the literature. We based these codes on both ASE studies and DNA methylation studies (supplementary table S3 and references therein, Supplementary Material online).

We generated imprinting codes at the SNP level after considering all genes affected by the polymorphism using VEP, Variant Effect Predictor (McLaren et al. 2016), and filtering for genes expressed at >0.1 TPM in a placental reference RNA-seq data set. The reference data were from a study with high quality RNA from Caesarean deliveries (Majewska et al. 2017). SNPs that mapped on the same strand to multiple genes with different imprinting codes were not included in our analyses. The proportion of paternal alleles was calculated for a subset of SNPs (35%) for which we had phasing information (N = 3,684, supplementary fig. S1, Supplementary Material online). The mean paternal allele bias for each group of SNPs was consistent with the literature in having the expected direction of imprinted expression, but it was also apparent that for many SNPs, allelic bias was only partial. The use of literature-derived imprinting codes enabled us to analyze MEGs, PEGs, and CEGs separately without the need for phased genotypes to experimentally determine parent of origin for each gene.

Statistical Analysis

We chose the exon as the unit of analysis as allelic bias along an exon should be constant. It is determined by the expression levels and allelic bias of the transcripts that share an exon (Deonovic et al. 2017). An exon was deemed to be imprinted if >90% of the reads for at least 50% of the SNPs within the exon mapped to one allele. These thresholds did not distort imprinting calls for exons with low numbers of SNPs or high expression (see supplementary analysis, Supplementary Material online) and accommodated outliers and technical sources of nonconcordance across an exon. Our dichotomous indicator of allele-specific expression is referred to as “ASE” below.

We used binomial GLMMs to produce a uniform set of person-level and gene-level imprinting rate estimates (Gelman and Hill 2006). These models yielded variance parameters that were used to assess the overall contribution of genes, exons, placenta, and cotyledon to imprinting rate variation. The imprinting status could be determined in each sample for only 15% of the exons in the targeting region. An advantage of the GLMM framework is that we were able to estimate these population-level variance parameters with incomplete data. A similar approach was used recently to analyze ASE from human brain RNA (Gulyas-Kovacs et al. 2018). The GLMMs included random effects for placenta and for cotyledon (person-level variables), and for genes and exons (genome-level variables). The random cotyledon effects were nested within the random placenta effects, and the random exon effects were nested in the random gene effects. The person-level and genome-level random effects were crossed with respect to each other. We also obtained best linear unbiased predictions (“BLUPs”) for all four random effects. In the final data set, we pooled data generated from multiple libraries (supplementary table S2, Supplementary Material online) taking advantage of the fact that ASE analysis compares read counts only within samples and thus is less affected by technical parameters than in studies of relative expression (Castel et al. 2015). Moreover, we included library type as a dummy variable in the GLMMs.

Maternal Contamination

The placenta is unique as it is composed of cells with different genotypes. Compared with mice, the maternal and fetal compartments in human placentas are less intermixed but it is still challenging to dissect a fetal sample that is completely free of maternal tissue (Benirschke et al. 2012). Maternal contamination will affect ASE in the fetal compartment of genes that are highly expressed in maternal decidua by increasing maternal reads. Thus, maternal contamination will overestimate LOI in PEGs and to a lesser degree underestimate it in MEGs. However, by quantifying nonfetal allele reads at SNPs that are homozygous in the fetus, a gene-level measure of contamination can be obtained directly from RNA-seq data (Hamada et al. 2016). Our analyses discarded genes contaminated above a threshold on a gene by sample basis (see Materials and Methods). All statistical models include a contamination covariate to control for any residual contamination (see Materials and Methods). We also used pyrosequencing and RT-qPCR to identify contaminated samples.

Results

Variable LOI for Many Genes

We quantified the extent of imprinting by measuring the fraction of exons that were silenced according to the definition given above, which permits up to 10% reads from the silenced allele and eliminates outliers. Our data have 5,478 data points (exons by sample) of which 39% had allelic bias of ≤90% indicating LOI. Here, we use the term LOI to indicate deviation of allelic bias from mono-allelic expression and imprinting rate is the dependent variable in the GLMMs. Because we did not attempt to measure imprints, LOI for a given gene does not necessarily correlate with the loss of an imprinting mark in the vicinity of the gene. When grouped by parental expression bias, 25% of the 2,386 PEG-exons had LOI, 43% of the 1,139 MEG-exons had LOI, and 55% of the 1,953 CEG-exons had LOI. A gene by placenta summary of the results is presented in supplementary figure S2, Supplementary Material online. This overview shows that LOI was widespread in maternally and paternally silenced genes.

The results from the GLMM analysis stratified by parental expression bias are presented in table 1.

Table 1.

Stratified Models.

Preferentially Expressed Allele	Paternal (PEG)				Maternal (MEG)				Complex (CEG)
Fixed Effects	Estimate	SE	Stand. Coeff.	P	Estimate	SE	Stand. Coeff.	P	Estimate	SE	Stand. Coeff.	P
Birth length cen. (cm)	−0.072	0.033	−0.138	0.031	−0.086	0.045	−0.161	0.057	−0.109	0.032	−0.206	0.001
Placental weight cen. (kg)	−4.55	0.738	−0.374	<0.001	−1.375	1.059	−0.109	0.194	−2.233	0.763	−0.182	0.003
Boya	0.003	0.123		0.981	−0.075	0.168		0.655	0.009	0.126		0.944
Later bornb	0.281	0.2		0.159	0.542	0.263		0.039	0.57	0.178		0.001
Gene class=lncRNAc	2.68	0.24	0.974	<0.001	3.404	0.319	1.353	<0.001	0.886	0.2	0.275	<0.001
Random Effects	SD	2.5% Bound		97.5% Bound	SD	2.5% Bound		97.5% Bound	SD	2.5% Bound		97.5% Bound
Placenta	0.747	0.614		0.909	0.061	0.050		0.075	0.043	0.035		0.052
Cotyledon	0.053	0.045		0.062	0.681	0.583		0.795	1.175	1.003		1.376
Gene	1.272	1.017		1.590	1.225	0.906		1.655	0.049	0.037		0.064
Exon	0.634	0.572		0.702	0.435	0.376		0.503	1.489	1.328		1.669
	N				N				N
Placentas	52				50				49
Cotyledon	86				83				80
Imprinting calls	2,324				1,114				1,915
Genes	40				22				29
Exons	192				95				153

Note.—Dependent variable is imprinting rate. Technical variables included as fixed effects but not shown: Library type (4 levels), Sequence batch (6 levels), RIN, Percent mitochondrial sequences, maternal contamination.

Ref. Category: girl.

Ref. Category: first born.

Ref. Category: protein coding.

As revealed by the random effects, both PEGs and MEGs showed strong gene-to-gene variation in ASE that exceeded the variation among exons in the same gene (fig. 1). Allelic bias is determined by the relative activity of the promoter that controls transcription of the gene on the maternal and paternal alleles. The finding of greater variation in ASE among genes than among exons is consistent with a single promoter driving the expression of most exons in a gene. Nonetheless, the variation at the exon level was substantial, suggesting that transcripts with unique combinations of exons from the same gene have distinct allelic biases due to the activity of their individual promoters. In contrast to MEGs and PEGs, in CEGs the variance component for exons exceeded that for genes. In the literature (supplementary table S3, Supplementary Material online), transcripts for the same CEG tend to have different expression biases, or studies disagree with respect to a CEG’s parental bias, or the parental expression bias was not specified. In our GLMM analysis, this complexity resulted in high variability for the exon random effect.

. 1.

Variance components for genes and exons in the GLMM models (table 1) by gene group. Error bars correspond to the 2.5% lower and 97.5% upper bounds.

The best linear unbiased prediction values (BLUPs) showed that both PEGs and MEGs lost imprinting (fig. 2 and supplementary table S4, Supplementary Material online). A BLUP is the log odds ratio that a gene is imprinted after adjusting for the fixed effects in the model; strongly positive BLUPs indicate that the gene has retained imprinting. Among the PEGs, some genes (RHOBTB3 and SNU13) showed frequent LOI while others (KCNQ1OT1, IGF2-AS, and PEG3) were rarely expressed from the maternal allele. Among the MEGs, TFPI2 stood out for its low imprinting rate and its high expression in placentas. Conversely, MEG3, MEG9, and H19 rarely lost imprinting.

. 2.

Gene-level imprinting scores derived from gene-class fixed effect and gene-level BLUPs. Higher values indicate a higher propensity for the gene to be imprinted (mono-allelic expression). Genes belonging to the gene class lncRNA are labeled by asterisk. Numbers to right indicate the number of exons with data for each gene in all samples. Crude values refer to average of all imprinting rates in all samples without adjustment for covariates in the GLMMs. Error bars are BLUP SDs.

To test whether the lncRNA genes as a class had lower LOI, we grouped genes into protein coding, lncRNA, and ncRNA and then added these indicators as fixed effects in the GLMMs (table 2). The lncRNA gene class was much more likely to be imprinted relative to protein coding genes (log odds ratio 2.874 P < 0.001).

Table 2.

Nonstratified Model (PEG+MEG).

Fixed Effects	Estimate	SE	Stand. Coeff.	P
Birth length cen. (cm)	−0.085	0.027	−0.161	0.001
Placental weight cen. (kg)	−3.457	0.604	−0.281	<0.001
Boya	−0.072	0.099		0.467
Later bornb	0.31	0.158		0.05
lncRNAc	2.874	0.191	1.078	<0.001
ncRNAc	−0.793	0.623		0.203
MEGd	−0.705	0.102	−0.33	<0.001
				N
Placentas				52
Cotyledon				86
Imprinting Calls				3,438
Genes				62
Exons				287

Note.—Dependent variable is imprinting rate. Technical variables included as fixed effects (not shown): library type (4 levels), sequence batch (6 levels), RIN, percent mitochondrial sequences, maternal contamination. Random effects similar to table 1 were included in the model (not shown).

Reference category: Girl.

Reference category: First born.

Reference category: Protein coding.

Reference category: PEG.

In sum, we found LOI in both PEGs and MEGs. Some genes were more prone to LOI than others and the variability in imprinting rates was greater at the gene level than at the exon level. In contrast to protein coding genes, lncRNA genes usually retained full imprinting.

Greater LOI in Longer Newborns

As imprinted genes are associated with growth phenotypes (Adkins et al. 2010; St-Pierre et al. 2012; Hoyo et al. 2014; Kappil et al. 2015; Moore et al. 2015), LOI has been suggested to contribute to phenotypic plasticity but experimental evidence for this concept has been lacking (Koukoura et al. 2012). The placentas that we analyzed were from an undernourished population with about half of the F1 mothers having experienced early childhood stunting (Strassmann 2011) and their F2 offspring having an average birth weight of 2.6 kg ±0.42, indicative of reduced fetal growth. Thus, we had the opportunity to investigate the hypothesis that LOI is functional and mediates phenotypic effects. We focused on birth length as it is a parameter that presages later stunting (MAL-ED Network Investigators 2017). The evidence for a negative impact on adult health and for transgenerational transmission is stronger for birth length than for birth weight (Victora et al. 2008; Martorell and Zongrone 2012). Birth length and birth weight were correlated (r = 0.56) and a larger sample size would be required to add both as covariates in one GLMM model. We also included placental weight in the models as some imprinted genes differentially affect fetal growth and placental weight (Moore et al. 2015).

To assess variation in imprinting rates across subjects, we predicted person-specific imprinting rates by summing the random effect BLUPs for placentas and cotyledons within an individual (supplementary table S5, Supplementary Material online). By doing so, we are able to include the repeatedly measured individuals due to sequences from multiple cotyledons of the same placenta (n = 18) or due to multiple sequences of the same cotyledon (n = 8). The difference in predicted imprinting odds ratio between extreme samples (most imprinted relative to least imprinted) were large (PEGs 3.1, MEGs 3.7, CEGs 12.3), indicating substantial variation in the imprinting rates among individuals. To find out if this variation was associated with birth length and placental weight, we included these parameters as covariates in the GLMMs. In the stratified analyses (table 1), all imprinting classes (PEGs, MEGs, and CEGs) had a negative estimate for birth length and PEGs and CEGs had an even stronger negative estimate for placental weight. These results were reproduced in the nonstratified analysis (PEGs+MEGs) where birth length and placental weight had negative estimates of −0.085 and −3.457, both with P ≤ 0.001 (table 2). These negative estimates indicate that newborns of increased birth length had less imprinting in their placentas and heavier placentas also had less imprinting.

LOI Did Not Correlate with Sex

In the GLMMs, the effect size for sex was small relative to its SEs; thus, we found no evidence that sex was associated with imprinting rates in any of the gene groups. This lack of association was not due to limitations of power as 42% of the placentas were from boys. In contrast, later borns did have increased allelic bias in MEGs and CEGs suggesting that parity modulates LOI. However, this effect was dependent on TFPI2 as deleting this MEG in the course of the sensitivity analysis resulted in a nonsignificant estimate for later borns (P = 0.39). Because only 15% of the placentas in this study were from later borns, reproducing these results in a larger data set is desirable. Deleting TFPI2 did not alter any of the other significant associations.

Paternally Imprinted Genes Had More LOI than Maternally Imprinted Genes

The stratified analyses revealed that both maternally and paternally imprinted genes lost imprinting. To find out if LOI was more prevalent for PEGs or MEGs, we added the fixed effect MEG to the unstratified GLMM. MEGs had a negative fixed effect (−0.705) with P < 0.001, indicating that they more frequently lost imprinting. To ensure that the preferential LOI of MEGs was not driven by a few genes or samples, we performed a sensitivity analysis by sequentially removing the four genes or the four samples with the most data. Removing genes reduced estimates for MEGs in the unstratified model by <10%; removing samples also did not change the estimates, the overall trend, or the significance, confirming that MEGs in general have more expression from the silenced allele.

LOI of MEGs and PEGs Was Correlated

As we observed significant LOI in genes with imprints from either parent, we were interested in finding out whether LOI of PEGs and MEGs was correlated in the same placenta. We calculated a correlation coefficient of 0.47 (95% confidence interval: 0.22–0.66). This degree of correlation suggests that when LOI occurs it affects multiple genes irrespective of whether their expression bias is maternal or paternal. Given that the silenced alleles of MEGs and PEGs reside on different chromosomes, we infer that the correlated modulation of allelic bias involves a trans-acting factor and is not regulated in cis.

Discussion

We analyzed placental specimens from a longitudinal and multigenerational cohort study to systematically characterize variation in LOI across genes and individuals. Since LOI often occurs in cancers, it has been associated with disease (Jelinic and Shaw 2007). Whether LOI has any functional consequences in healthy tissue has yet to be established. The phenotypes of mice with deleted imprinted genes and current understanding of the evolutionary forces that led to imprinting (Peters 2014) suggest that LOI may fine-tune fetal growth in placentas and additional physiological processes in other tissues (Tucci et al. 2019). Thus, we tested the hypothesis that modulation of imprinting is an epigenetic mechanism for regulating offspring growth, as measured by birth length and placental weight. We show that newborns of increased birth length and heavier placentas had more LOI. Our results contradict the alternative view that LOI is a stochastic process without phenotypic consequences.

LOI in the Placenta Is Common

Over one-third (39%) of the exons in the data set had some expression from the silenced allele revealing that LOI is common in many imprinted genes. Intriguingly, imprinted genes that code for lncRNAs rarely lost imprinting (fig. 2). These genes, where known, function mainly as epigenetic regulators of multiple targets (Barlow 2011). The maintenance of mono-allelic expression may indicate that gene dosage of master regulators is tightly controlled and LOI of lncRNA is incompatible with a viable fetus. The targets of lncRNAs that are also imprinted showed a wider range of allelic bias (supplementary table S4, Supplementary Material online). For example, KCNQ1OT1 is a lncRNA that had a BLUP of +3.1 (implying mono-allelic expression) and it controls the imprinted protein coding genes KCNQ1, CDKN1C, and PHLDA2 that had BLUP values of +0.04, +0.52, and +1.00, respectively. The variability of LOI in the targets suggests that suppression by lncRNA is incomplete and that additional factors influence allelic bias.

Frequent deviation from complete allelic bias was also found in the multiple tissues collected as part of the GTEx project, in which placentas were not included (Babak et al. 2015; Baran et al. 2015), and in a study of the human cortex (Gulyas-Kovacs et al. 2018). Allelic bias is also subject to coordinated spatiotemporal regulation as was shown in mouse brains (Perez et al. 2015). LOI was not due to maternal contamination as we excluded genes from placentas contaminated with maternal RNA and genes with unknown contamination levels (see Materials and Methods). Clustering by ASE Spearman correlation showed that samples from two cotyledons of the same placenta clustered together (rs = 0.50) while samples from different placentas did not (rs = 0.18) (supplementary fig. S3, Supplementary Material online). Thus, LOI was consistent between two sampling sites on the same placenta. Similarly, DMRs from different biopsies showed higher correlation in methylation when they originated from the same placenta (Monteagudo-Sánchez et al. 2019). In sum, derepression of imprinted alleles is not limited to disease states but also occurs in many healthy tissues—including the placenta—as part of a normal developmental program.

Stronger LOI for MEGs than PEGs

We observed that MEGs lost imprinting more frequently than PEGs, suggesting that paternal imprints are removed more frequently than maternal ones. As MEGs are silenced by imprints on paternal alleles, LOI of these genes skews allelic bias in the placenta toward expression of growth inhibiting genes and thus toward the interests of the alleles derived from the mother.

The mechanisms that lead to LOI of imprinted genes are likely to involve loss of the imprint, not acquisition of new repressive marks on the preferentially expressed allele as most genes maintain the direction of their allelic expression bias during development into adulthood (Wilkins et al. 2016; Zink et al. 2018; Monk et al. 2019). Currently, the only imprint characterized during development is DNA methylation. DNA methylation is established during germ cell development but DMRs only emerge during preimplantation development as a result of their resistance to demethylation (Hanna and Kelsey 2014). The discovery of factors that protect from demethylation shortly after fertilization and that modify multiple DMRs implicate the preimplantation period as critical for shaping imprints (Payer et al. 2003; Mackin et al. 2018). Some of these factors that protect DMRs are primate specific (Takahashi et al. 2019). Furthermore, this period prior to zygotic gene activation is also longer in primates compared with rodents and involves the establishment of a species-specific chromatin conformation (Xia et al. 2019). The involvement of maternal transcripts in the reactivation of MEGs during this period had already been postulated by Wilkins and Haig (2002) through a kinship model that included parental trans-acting modifiers of imprinting.

Given these mechanistic insights and the association of LOI with birth length, we propose that LOI of multiple imprinted genes is a mechanism that adjusts offspring size to maternal resources, in part through trans-acting maternal modifiers that act prior to the maternal to zygotic transition. Both maternal and paternal imprints are being remodeled during this period such that paternal imprints that represent the interests of the paternal genome are being modified by the mother, who usurps control.

Correlated LOI of MEGs and PEGs

Some individuals had higher LOI in many genes as shown in the BLUP tables that display a wide range in the combined placenta+cotyledon BLUPs (supplementary table S5, Supplementary Material online). LOI in MEGs and PEGs was correlated, which implies a degree of coordinated regulation of LOI.

According to the kinship hypothesis PEGs and MEGs are functionally antagonistic, so coordinated LOI is an unexpected finding and would be predicted to reduce the phenotypic consequences of LOI. However, it is possible that the two gene groups differ in their phenotypic penetrance. In our data, the association of increased birth length and placental weight with LOI indicates that less LOI in the more numerous PEGs had a greater effect on the phenotype than more LOI in the fewer MEGs. Thus, the correlated LOI of MEGs and PEGs did not cancel out phenotypic effects. It is unclear if correlated LOI reflects a mechanistic constraint or a phenotypic constraint. In the latter case, LOI that is too asymmetrical may not result in a viable offspring. Greater clarity will require a better understanding of the quantitative relationship between LOI, expression, and placental function of imprinted genes.

While the prezygotic phase of development is likely to be the developmental stage during which LOI occurs, all the effector molecules that orchestrate correlated LOI are not known. However, multiple locus imprinting disturbances (MLID) provide evidence that mechanisms that affect multiple DMRs do exist (Sanchez-Delgado et al. 2016). Some of the genes that cause MLID are maternal effect genes, reinforcing the view that the developmental window for LOI is soon after fertilization (Mackay et al. 2008; Begemann et al. 2018).

Correlated LOI between PEGs and MEGs implies that some gene regulation operates through changes in repression of many imprinted genes. However, gene-specific levels of regulation are well established and are expected to degrade this correlation. For example, multiple IGF2 promoters have been identified that produce transcripts with varying tissue specificity (Monk et al. 2006). Tissue-specific promoters can even lead to change in the parental origin of the imprint as is well documented for GRB10 (Monk et al. 2009) and IGF2 (Monk et al. 2006; Baran et al. 2015).

Positive Association between LOI and Offspring Birth Length and Placental Weight

The majority of studies have used DNA methylation to assess associations between LOI and phenotypic outcomes. A recent meta-analysis of these studies concluded that there were no consistent associations in part due to methodological issues (Maddock et al. 2018). Studies focused on the methylation of a few selected imprinted genes in placentas also did not find consistent associations (McMinn, Wei, Schupf et al. 2006; Guo et al. 2008; Bourque et al. 2010; Koukoura et al. 2011), which is likely due to the incomplete understanding of the relationship between DNA methylation and expression as well as the possibility of other imprints contributing to the regulation of allelic bias (Inoue et al. 2017). Studies that have analyzed relative expression of some imprinted genes have found correlations consistent with the kinship hypothesis (Apostolidou et al. 2007; Guo et al. 2008; Lambertini et al. 2012), particularly from first trimester placentas (Moore et al. 2015). Relative RNA expression levels from term placentas can be dominated by the birthing process itself (Lee et al. 2010) or other transient exposures. Measuring changes in allelic bias, as we did here, allowed us to use term placentas to estimate ASE, an RNA property that is under epigenetic control and arguably less affected by short-term stimuli.

In PEGs and CEGs, the association between placental weight and imprinting rate was much stronger than that for birth length. Although the placenta is a discarded organ, placental weight may have many downstream effects and is known to be associated with adult health (Barker et al. 1990). Longer newborns had more biallelic expression in their placentas and placentas with more biallelic expression were heavier. Placentas of shorter newborns tended to express PEGs and CEGs from one allele, consistent with the proposal that resource scarcity produces a greater degree of genomic conflict and results in higher allelic bias. The associations between imprinting rate and birth phenotypes supports our hypothesis that modulation of LOI is functional and agrees with the hypothesis that PEGs are growth promoting, particularly in the placenta.

For MEGs a positive association between imprinting rate and placental weight may be expected as these genes are predicted to be growth inhibitory. However, the estimate was negative and not significant, thus arguing against a role for LOI of MEGs in reducing placental weight and birth length even in conditions where half of the mothers were stunted. In addition, our data are also not compatible with a compensatory role for LOI in which shorter fetuses upregulate the expression of growth promoting PEGs by activating the silenced allele.

The concept of an imprinted gene network has been proposed to explain coregulation of many genes that control embryonic growth in mice (Varrault et al. 2006). This concept was developed from experimental manipulation of apical transcription factors in the network and therefore is different from what we observe here, which is an allele-specific repressive mechanism that is thought to operate on chromatin-level properties impacting imprinted genes.

We collected the placentas from women in a prospective cohort study. Half of the newborns had low birth weight and many were probably small for gestational age though conception dates were unavailable. Nonetheless, most infants thrived and had no overt developmental delays at age 1 year. Thus, the variation in LOI that we documented was seen in phenotypically normal children growing up in an environment with food insecurity. Currently, there are insufficient data available to determine whether LOI correlated with placental weight or birth length in other populations (Guo et al. 2008; Heijmans et al. 2008; Diplas et al. 2009; Einstein et al. 2010; Tabano et al. 2010; Tobi et al. 2011).

Because LOI is derived from the ratio of expression of two alleles, it does not uniquely define the changes in expression. Studies in the mouse cerebellum have shown that LOI most frequently resulted from reduced gene expression, especially from the preferentially expressed, nonimprinted allele (Perez et al. 2015). For many other genes, however, reduced allelic bias was the result of increased expression of the silenced allele or changes in the expression of both alleles (Perez et al. 2015). Similarly, there was no straightforward correlation in human tissues between the relative expression of a gene and allelic bias when comparisons were made across tissues (Baran et al. 2015). We were unable to compare relative expression between samples because of the variable quality of the RNA. Thus, absent relative expression estimates for each allele, we cannot infer what change in expression of which allele(s) produced LOI. Further, the short read length of the RNA-seq technology limited our ability to determine LOI at the transcript level. Despite these limitations, LOI is indicative of expression changes affecting one or both alleles.

This study is the first genome-wide survey of LOI in known imprinted genes in human placentas. Other epigenetic studies from West Africa exploited the seasonal availability of food to identify metastable epialleles in blood (Dominguez-Salas et al. 2014). One allele that stands out for its high significance scores and its association with birth season overlaps the transcription start site of the VTRNA2-1 gene. Remarkably, the methylation status of its DMR is stable for more than a decade after birth (Silver et al. 2015). It has yet to be established whether hypomethylation of this DMR leads to LOI and we could not determine LOI of VTRNA2-1 in our placental samples because expression of this gene was too low. Interestingly, a reanalysis of these data concluded that children born during the rainy season when food was “low calorie nutritionally rich” were more likely to be imprinted at this DMR (Carpenter et al. 2018). This report for VTRNA2-1 is consistent with our transcriptome-wide data-linking allelic bias to birth phenotypes.

In sum, we took advantage of the almost completed discovery of imprinted genes in humans to systematically characterize the variation in LOI of 91 genes expressed in placentas. Variation in LOI was primarily at the gene as opposed to the exon level and some genes were less susceptible to LOI than others—in particular, lncRNA genes rarely deviated from mono-allelic expression. We show that LOI was common and nonstochastic and highly variable between individuals. Within a placenta, the expression bias of maternally and paternally imprinted genes was correlated. We propose that the regulatory mechanisms that mediate correlated LOI act during the preimplantation period and involve trans-acting maternal effect genes. The kinship hypothesis identified imprinted genes in the placenta as critical for growth. Building on this insight, we hypothesized that LOI provides the mother with a means to fine-tune fetal size. In support of this hypothesis, we show that birth length and placental weight were associated with allelic bias, making this the first comprehensive report of an association between LOI and a birth phenotype.

Materials and Methods

Study Population

The mothers for the placental collection were participants in an ongoing multigenerational prospective cohort study in Mali, West Africa (Strassmann 2011). The participants (F1, N = 1,698) were born in the years 1993–2000 in a cluster of nine rural villages belonging to the Dogon ethnic group. We took anthropometric measurements upon enrollment (1998–2000) and at yearly intervals through 2000. The cohort was also measured in 2004 and 2007 and annually from 2010 to 2019 when the participants were in puberty and young adulthood. The F2 (N = 707 as of September 1, 2019) are being measured every 3 months. The F0 generation was measured once (N = 1,217).

IRB

Informed consent or assent was obtained from participants depending on whether they were adults or children. IRB approval was obtained from the University of Michigan IRBMED (HUM00043670) and from La Faculté de Médecine de Pharmacie et d’Odontostomatologie (FMPOS) de Bamako in Mali (No. 2016/68/CD/FMPOS).

Sample Collection

All samples were collected at the local hospital. Its maternity clinic is staffed by midwives and supported by a physician. The midwives identified subjects with the help of a roster of all women enrolled in the study who had given their informed consent for the placental collections. The same midwives performed prenatal consultations and they were in most cases already familiar with the women prior to the day of parturition. To insure that the health of the baby and the mother were not compromised, the collections were always carried out by a different midwife from the one attending to the birth.

On the day of collection, formaldehyde (37%) was diluted 1:10 with phosphate buffered saline and stored on ice. The diluent was prepared with a PBS tablet (Sigma–Aldrich) and distilled water. In preparation for the birth, RNAlater aliquots (4.5 ml, Sigma–Aldrich) in 5 ml cryovials stored at −20 °C and Tempus Blood RNA tubes (Thermo Fisher Scientific) were transferred to the maternity clinic on ice. Five wash tubes with 15 ml PBS were precooled on ice as well.

The local best practice protocol for singleton births is to accelerate the delivery of the placenta with oxytocin (10 IU, IM) immediately after delivery of the baby. All specimen collections were done within 30 min of the recovery of the placenta. First, 3 ml of umbilical cord blood was collected from the umbilical arteries with a syringe equipped with a 21 Gauge needle, added to the Tempus tube, and put on ice. Second, the placenta was freed from the fetal membranes and placed with the umbilical cord facing down on the dissection tray. In a circle encompassing a central area that is three quarters of the placental surface, two well-formed cotyledons were identified. Each cotyledon was cut in half and ∼1 g of tissue was dissected from the fetal compartment in the interior. The surface of each cotyledon was removed with a shaving motion of the scalpel to yield the maternal decidua (not analyzed here). The tissue samples were minced finely, washed in cold PBS, and transferred to 4.5 ml of RNAlater. Samples were incubated in RNAlater at 4 °C for 36 h before excess liquid was drained and the samples transferred to a −20 °C solar freezer. Fixed specimens of the fetal compartment were obtained from the opposite cut face at the center of each cotyledon. They were immediately placed in tissue cages and submerged in 3.7% formaldehyde for 36 h, followed by a 30 min wash with 70% ethanol and storage at −20 °C. Samples were stored up to 10 months at −20 °C prior to shipment on dry ice via World Courier to a laboratory freezer (−80 °C) at the University of Michigan.

Nucleic Acid Purification

RNA

The TRIzol protocol was used with the following modifications (Chomczynski 1993). Tissue homogenization was achieved in two steps. First, ∼50 mg of thawed placental tissue was crushed frozen in a prechilled (LN2) cell crusher (cellcrusher.com) with 30 hammer strokes. Second, the tissue powder was transferred into 1 ml TRIzol and homogenized with 3×30 s bursts of a TissueRuptor (Qiagen). Phase separation was achieved with 100 μl 1-bromo-3-chloropropane (Sigma–Aldrich). The RNA was immediately further purified with a Zymo-Spin Column (R1016, Zymo Research) and digested with DNaseI for 30 min at 37 °C (AM1907, Ambion).

DNA

Genomic fetal DNA was isolated from umbilical cord tissue with QIAamp DNA Mini Kit following the manufacturer’s instructions (Qiagen).

DNA from TRIzol preparations

To measure maternal contamination by pyrosequencing, DNA was isolated from the same tissue sample as the RNA by recovering DNA from the organic phase following a standard protocol (Stephenson et al. 2016).

Pyrosequencing

Validation of RNA-seq allelic bias was performed on select SNPs using pyrosequencing. cDNA synthesis RT was performed immediately after DNaseI digestion of RNA with the ProtoScriptII First Strand cDNA Synthesis Kit (E6560, New England Biolabs) and random hexamer primers. Pyrosequencing primers were designed with Qiagen PyroMark Assay Design 2.0 software and amplicons were generated with PyroMark PCR Kit (978705, Qiagen) and sequenced using a PyroMark Q96 MD workstation. 30 SNPs were pyrosequenced in 22 cotyledons and the results are shown in supplementary table S6, Supplementary Material online. Pearson correlation between LOI determined by RNA-seq and pyrosequencing was 0.59 (P = 9.2×10−9, n = 90) for the SNPs that were included in the statistical model.

Library Preparation

Libraries were prepared by the University of Michigan DNA Sequencing Core following the manufacturer’s instructions. All samples were spiked with ERCC control RNA (mix 1 or mix 2) according to the manufacturer’s instructions (4456739, ThermoFisher Scientific). In most cases the library kit used was SMARTer Stranded Total-Seq Kit Pico Input Mammalian (635006, Clontech) using 10 ng of input RNA. Samples with RIN >6 were also sequenced with the High Mammalian version of the same kit (634873, Clontech) using 1 μg input RNA. Eleven sequence libraries were prepared with NEBNext Ultra Directional (E7760, NEB) using 1 μg input RNA, and two with Illumina TruSeq RNA Access (RS301, Illumina) using 1 μg input RNA. See supplementary table S2, Supplementary Material online, for summary statistics.

RNA Sequencing

Paired-end sequencing at 150 bp read length was performed for the bulk of the samples on an Illumina HiSeq 4000 platform at a depth of one lane per sample. Read length was 75 bp for some lower RIN samples and four samples were multiplexed on one lane for NEB and TruSeq libraries.

DNA Sequencing

Targeted DNA sequencing was employed for genotyping using the Roche-NimbleGen design and SeqCap EZ Choice. The targeted region included all exons extended at each end by 50 bp of imprinted genes and 1 kb sequence at the 5′ end (0.99 Mb). Genes are listed in supplementary table S3, Supplementary Material online. In addition, custom selected genes were added bringing the total targeted region to 1.42 Mb. KAPA HyperPrep libraries (KK8504, Roche) and enrichment were performed according to the manufacturer’s instructions and samples multiplexed no more than 72-fold and paired-end sequenced on one HiSeq 4000 lane at 125 bp read length.

DNA-Sequence Analysis

Illumina adapter contamination and read ends with base quality <20 were removed using Trimmomatic v0.36 (Bolger et al. 2014). Reads shorter than 75 nt after trimming were discarded. Trimmed reads were aligned to hg19 reference genome using BWA v0.7.15 (Li and Durbin 2009). Duplicate reads were removed using MarkDuplicated from Picard v2.7.1 (http://broadinsitute.github.io/picard). GATK v3.7 with its hg19 resource bundle (Van der Auwera et al. 2013) was used to perform base quality recalibration by BaseRecalibrator and variant calling by HaplotypeCaller in its joint-calling mode. Resulting variants underwent GATK-recommended hard-filtering for SNPs and INDELs separately. Variant calling and filtering were restricted to the targeted region padded by 500 bp. Furthermore, we applied a series of filters in order to remove less-confident genotypes that included the following: 1) sites with genotyping quality <20 or total read depth <20; 2) sites that fell in regions with 75mer-alignability score <1 using the ENCODE mappability track (Derrien et al. 2012; Rosenbloom et al. 2012); 3) sites in the DAC (ENCODE Data Analysis Consortium) Blacklisted Regions or the Duke Excluded Regions (ENDCODE Project Consortium 2012); 4) variants with known alternate allele mapping bias identified in this study (Panousis et al. 2014; Castel et al. 2015); 5) SNPs that were called as INDELs in another sample; 6) less-frequent alternate alleles at the SNPs that had more than one alternate allele in the batch; 7) heterozygous SNPs whose reference allele frequency was <0.2 or >0.8 (Nielsen et al. 2011); 8) homozygous SNPs whose reference allele frequency was >0.05; 9) homozygous reference sites whose reference allele frequency was <0.95; and 10) SNPs where >5% of reads supported an allele that was neither reference nor alternate. In order to apply the above filters, allele-specific read depth at SNPs was recalculated using ASEReadCounter in GATK. PhaseByTransmission in GATK was used to phase the variants in F2 samples whose parents were both genotyped. The phased variants were filtered by requiring the transmission probability score to be no lower than 20, and then combined with the variants phased by HaplotypeCaller. Samples with excessive Mendelian violations indicative of nonpaternity or tube error were excluded from the final phasing results. The coordinates of the final variants were lifted over to hg38 using LiftoverVcf from Picard.

RNA-Sequence Analysis

Illumina adapter contamination and read ends with base quality <20 were removed using Trimmomatic v0.36. The first six nucleotides of the reads were trimmed off as suggested by the SMARTer kit manual. Reads shorter than 36 nt after trimming were discarded. HISAT2 v2.1.0 (Kim et al. 2015) was used to first build a new reference for each individual to incorporate the genomic variants identified from the corresponding DNA sample, and then to align the paired trimmed reads onto this reference with GENCODE GTF (release 27) (Harrow et al. 2012) as the known splice sites. Alignments were filtered using WASP v0.2 (van de Geijn et al. 2015) to reduce allelic bias. Then properly paired alignments with the highest mapping quality were selected as confident alignments and used for downstream analyses. StringTie v1.3.4d (Pertea et al. 2016) was used to quantify the relative expression at the gene level. Alignments were split into sense-strand and antisense-strand alignments, and ASEReadCounter in GATK was used to calculate allele-specific RNA read depth in both strands, at each heterozygous SNP in the proband individual identified from the paired DNA sample. Parameter “–countOverlapReadsType” was set to “COUNT_FRAGMENTS_REQUIRE_SAME_BASE” to count only once in a read pair and require the two reads to be consistent at the examined base. Parameters “–minMappingQuality” and “–minBaseQuality” were set to 40 and 15, respectively, to exclude ambiguous bases. SNPs were annotated with the coordinates of the exons to which they mapped and with the names of the genes to which the exons belonged. Overlapping exons in the same gene were merged into one interval.

Assessment of Maternal Contamination

Maternal contamination in RNA-seq samples was assessed by examining the genotype at positions that are homozygous in the DNA but appeared to be heterozygous in the RNA (Hamada et al. 2016). All genomic locations in the targeted region with at least one heterozygous SNP in the population were examined. The fraction of non-Ref reads was calculated for these locations when the genotype of the fetus derived from umbilical cord tissue was homozygous Ref (0/0). A gene by sample contamination variable was obtained by averaging over all these locations that map to the same gene in a given sample. For a subset of samples (60%) for which the genotype of the mother was known, we determined the contamination similarly except that only locations were considered where the mother was heterozygous yielding a direct measure of the fraction of contaminating alleles. This revealed that 89% of the gene-sample pairs in this subset had a nonreference allele frequency <0.05. The full data set was then filtered on the 89th percentile of the contamination variable that was based only on the genotype of the fetus. In addition, we excluded from the analysis any gene by sample pairs for which there were no data to determine contamination.

To uncover any residual effects of maternal contamination after eliminating gene sample pairs, we included a sample-level contamination variable in the GLMMs as a fixed effect. This variable is the average frequency of non-Alt reads at homozygous Alt. SNPs (1/1).

The SNP-level sequencing data were converted to exon-level data and analyzed using the GLMM framework implemented in Python (https://github.com/statsmodels/statsmodels/blob/master/statsmodels/genmod/bayes_mixed_glm.py). The dependent variable is imprinting rate, which is determined from the imprinting status of each exon. Imprinting status is a dichotomous variable determined from all informative SNPs that map to the same exon. An exon was deemed to be imprinted if >90% of the reads for at least 50% of the SNPs within the exon mapped to one allele.

The following technical variables were included in the model as fixed effects: Library type: SMARTer Pico = Ref., SMARTER High Mammalian, NEBNext, Illumina TruSeq RNA Access; Sequence batch: 6 batches; GeneClass: Protein coding = Ref., lncRNA, ncRNA; RIN; Percent mitochondrial sequences; Maternal contamination.

Click here for additional data file.