Mapping the Cord Blood Transcriptome of Pregnancies Affected by Early Maternal Anemia to Identify Signatures of Fetal Programming.

CONTEXT: Anemia during early pregnancy (EP) is common in developing countries and is associated with adverse health consequences for both mothers and children. Offspring of women with EP anemia often have low birth weight, which increases risk for cardiometabolic diseases, including type 2 diabetes (T2D), later in life.
OBJECTIVE: We aimed to elucidate mechanisms underlying developmental programming of adult cardiometabolic disease, including epigenetic and transcriptional alterations potentially detectable in umbilical cord blood (UCB) at time of birth.
METHODS: We leveraged global transcriptome- and accompanying epigenome-wide changes in 48 UCB from newborns of EP anemic Tanzanian mothers and 50 controls to identify differentially expressed genes (DEGs) in UCB exposed to maternal EP anemia. DEGs were assessed for association with neonatal anthropometry and cord insulin levels. These genes were further studied in expression data from human fetal pancreas and adult islets to understand their role in beta-cell development and/or function.
RESULTS: The expression of 137 genes was altered in UCB of newborns exposed to maternal EP anemia. These putative signatures of fetal programming, which included the birth weight locus LCORL, were potentially mediated by epigenetic changes in 27 genes and associated with neonatal anthropometry. Among the DEGs were P2RX7, PIK3C2B, and NUMBL, which potentially influence beta-cell development. Insulin levels were lower in EP anemia-exposed UCB, supporting the notion of developmental programming of pancreatic beta-cell dysfunction and subsequently increased risk of T2D in offspring of mothers with EP anemia.
CONCLUSIONS: Our data provide proof-of-concept on distinct transcriptional and epigenetic changes detectable in UCB from newborns exposed to maternal EP anemia.

Anemia is a major health problem globally, especially among pregnant women. About 42% of pregnant women worldwide suffer from anemia defined as hemoglobin (Hb) < 11 g/dL in the first and third trimesters and < 10.5 g/dL in the second trimester (1). The prevalence of anemia during early pregnancy (EP) is about 18% in developed countries, but 35% to 75% in developing countries (2). Sub-Saharan Africa has a high incidence of anemia and in Tanzania, the prevalence was reported to be 57% in 2015-2016 (3). Anemia can be caused by several factors including malnutrition, infections such as malaria, and hereditary disorders (4-6). In Tanzania, malaria together with iron deficiency are common causes of anemia during pregnancy (7-9).

Anemia during pregnancy can have severe effects on the health of the offspring. Moderate to severe anemia (Hb < 8-9 g/dL) has been associated with decreased birth weight (BW) (10, 11), and first and second trimester anemia may be especially detrimental (12-14). Barker and Hales hypothesized that insults during fetal life causing decreased BW could predispose an individual to conditions like diabetes and cardiovascular disease (CVD) later in life (15-18). How and when anemia affects fetal growth is not well established, but poor growth prior to the third trimester of pregnancy could represent a particularly vulnerable time-window of pregnancy for fetal programming predisposing to diabetes and CVD (19). In Tanzania, anemia and infectious diseases remain high in children, while diabetes and CVD increase rapidly in older generations (20, 21), facilitating a need to explore the role of maternal anemia on the epidemic of diabetes and CVD.

The impact of anemia during EP on the risk of chronic diseases in the offspring may potentially be mediated by epigenetic mechanisms, for example, DNA methylation, which can influence gene expression in the offspring, here represented by umbilical cord blood (UCB). Therefore, examination of gene expression alterations in the UCB from offspring of mothers who had anemia early in pregnancy could provide insights on how this sensitive fetal period may predispose the offspring to chronic disease later in life. To address these questions, we performed RNA sequencing of UCB from offspring of Tanzanian mothers with EP anemia from the FOETALforNCD study (22). We furthermore studied the accompanying changes in epigenetic patterns among selected genes whose expression levels were altered by early maternal anemia. Subsequently, we assessed if the expression levels of these genes were associated with short and/or long-term outcomes in the offspring. We studied the expression levels of the same genes in fetal and adult liver as well as pancreatic samples from another study (Swedish, GTEx) to determine if they had any potential roles in development and/or function of these organs. We compared insulin levels in UCB from EP anemic mothers to controls to identify early signs of impaired insulin secretion or action among EP anemic offspring. We then tested if any of the genes whose expression was altered in UCB of early-anemic mothers could have an impact on beta-cell development/function, which could then explain the lower cord insulin levels and higher type 2 diabetes (T2D) risk in later life.

Methods

Study Participants

The FOETALforNCD cohort profile has been described in detail previously (22). In brief, the field study was conducted in a rural region of Northeast Tanzania, in Korogwe District and adjacent villages in Handeni District, Tanga Region, from July 2014 to December 2016. The central field site was the maternity ward and the reproductive and child health clinic (RCH clinic) at Korogwe District Hospital and mobile outreach clinics were set up in 48 villages. Participants were recruited prior to conception or in early pregnancy and followed throughout pregnancy until birth and 1 month postnatally. In total, 1415 nonpregnant women participated in the preconception study, 538 in the pregnancy study, and 427 in the birth cohort study. Data collection included: maternal blood, screening for noncommunicable diseases and malaria, ultrasound in each trimester, neonatal anthropometry at birth and at 1 month of age, and UCB (22). Briefly, UCB samples from 366 participants were collected in EDTA, serum, and PAXgene tubes as described previously (22). The study was approved by regional ethical committees and all participants gave written informed consent.

Study Design and Maternal Diagnosis

Maternal hemoglobin (Hb) was measured at the time of inclusion, all antenatal visits, and at delivery. Inclusion criteria for the main cohort in the presented analyses were: availability of both UCB and UCB collected in PAXgene tubes (for RNAseq), malaria-negative throughout pregnancy, BW measured within 24 hours of delivery, gestational age of > 22 weeks at delivery (≥ 155 days), HIV negative (unknown HIV status excluded), a live birth, not having any hypertensive disorders, not having high Hb (>14 g/dL) at any point during pregnancy, and only singleton births (ie, no twins or triplets).

Anemia cases

Anemic cases comprised 50 women with early anemia defined as Hb ≤ 10.9 g/dL at the earliest time point during the first and second trimesters. We selected 48 cases based on those who had the lowest Hb levels at the lowest gestational age and who were consistently anemic during early pregnancy, and reached 50 cases when the cutoff at a gestational age of 21 weeks + 5 days was set (Supplementary Table S1 (23)).

Controls

The 50 controls were selected on the condition that they (1) did not have anemia up to 21 weeks + 4 days of gestational age (GA) during all the Hb measurements, (2) Hb was within normal range (ie, Hb not > 15, since Hb levels up to 15 are within normal range for women; all mothers had Hb < 14, with the exception of 1 mother who had a Hb = 14.4 only at delivery and was included in the study) (3) did not show any drop of an empirically selected > 2.3 g/dL from the lowest Hb in the GA interval (<21 weeks + 4 days) and until delivery, (4) did not, at any time point in pregnancy, have an Hb < 10.0 g/dL. Based on these criteria, 50 women with the highest Hb levels were selected as controls (Supplementary Table S1 (23)).

RNA Isolation and DNA Isolation and Quantification

DNA was extracted using the QIAamp DNA Blood Midi Kit (Qiagen, catalog # 51185, Sollentuna). DNA was quantified using picogreen assay (Promega). RNA was extracted from cord blood stored in PAXgene tubes using the PAXgene Blood RNA Kit (Qiagen, catalog# 762174, Sollentuna). RNA quantity was assessed using the Agilent 4200 Tape station system and quantified using the Qubit assays (Thermofisher Scientific).

RNA Sequencing

Library preparation, sequencing, and alignment

A quantity of 1 μg of RNA with RNA integrity number (RIN) ≥ 7 was used for library preparation using the TruSeq Stranded Total RNA Library Prep Globin (Illumina). RNA sequencing was performed using the 75-bp paired-end protocol on a Nextseq 500 platform (Illumina). Each individual transcriptome yielded, on average, 24.3 ± 10 (mean ± SD) million paired-end reads mapped to the human genome (NCBI build 38). Paired-end 101-bp length output reads were aligned to the human reference genome (hg19) using STAR (2.4.1a) (24) with the following parameters:

STAR –genomeDir ../Human/STAR_Genomes/hg38_Gencode22 –sjdbGTFfile genes.gtf –readFilesIn file.fastq.gz –out FilterType BySJout –runThreadN 12 –outFilterMultimapNmax 20 –alignSJDBoverhangMin 1 –sjdbOverhang 99 –outFilterMismatchNoverLmax 0.04 –alignIntronMin 20 –alignIntronMax 1000000 –outSAMstrandFieldintronMotif –outReadsUnmappedFastx –readFilesCommandzcat

The annotated RefSeq GTF transcript and fasta genome files were downloaded from Ensembl (https://www.ensembl.org/index.html). Gene expression was measured as the normalized sum of expression of all exons. Exons were defined as nonoverlapping unique exonic units, as described previously (25). The dexseq_count python script (https://www.huber.embl.de/pub/DEXSeq/) was used by counting uniquely mapped reads in each exon. The htseq_count python script (http://www-huber.embl.de/users/anders/HTSeq/doc/counting.html) was used for counting uniquely mapped reads for each gene (26).

Preprocessing and statistical analyses

Gene and exon expression matrices were compiled after concatenating the gene and exon level counts, respectively, from all individuals and converting them to a list-based data object called a DGEList using edgeR (27). Expression was normalized for library sizes and expressed in terms of counts per million (cpm). Very low counts across all libraries provide very little information on differences in gene expression between different phenotypes; therefore, such genes were filtered out: the criteria of at least 5 cpm in at least 10% of the cases was retained. The quantile-adjusted conditional maximum likelihood (qCML) method was used to calculate the likelihood by conditioning on the total counts for each tag and using pseudo counts after adjusting for library sizes. The exact test based on the qCML methods was used to compute differences in gene expression in UCB between cases and controls adjusted for maternal age, gestational age at delivery, sex, batch, maternal body mass index (BMI) at inclusion, and UCB composition. UCB composition comprised cord Hb levels, T cells (CD8, CD4), B cells, large granular lymphocytes (NK cells), monocytes, and granulocytes derived from methylation profiles in UCB (28).

For assessment of correlation with phenotypes, 137 genes showing altered expression in UCB from newborns of early-anemic mothers were selected. Correlation of gene expression of the selected genes with maternal Hb levels in EP and maternal Hb levels at delivery was performed using the partial correlation test by applying the “ppcor” package in R (29). The correlations were adjusted for maternal age, gestational age at delivery, sex, batch, maternal BMI at inclusion, and UCB composition.

Association with placental weight, placental size, birth weight, waist-to-length ratio, and cord insulin and c-peptide levels was performed using Pearson correlation test after log normal transformation of the data in R or python.

Global DNA Methylation

Genome-wide DNA methylation analyses were performed on DNA extracted from buffy coats on the Illumina Infinium EPIC Bead Chip (WG-317-1001) with Infinium assay using the standard Infinium HD Assay Methylation Protocol Guide (part number 15019519, Illumina). Samples were randomly distributed on the arrays. For these assays, 1 µg of DNA was bisulfite-treated according to protocol (EZ DNA Methylation Kit, Zymo Research), hybridized to the Illumina Infinium EPIC Bead Chips and scanned on Illumina iScan according to protocol. The Illumina Infinium EPIC Bead Chip contains > 850 000 probes with 99% coverage of RefSeq genes with the capacity for 12 samples per chip.

The Genome Studio methylation module software of the GenomeStudio Genome Browser (NCBI build 38) was used to calculate the raw methylation score for each DNA methylation site, which is represented as the methylation β value. The β values were calculated as:

β = intensity of the methylated allele [M] ÷ (Intensity of the unmethylated allele [U] + intensity of the methylated allele [M] + 100)

All samples passed GenomeStudio quality control steps based on built-in control probes for staining, hybridization, extension, and specificity and displayed high quality bisulfite conversion efficiency with an intensity signal above 4000 (30). Probes with detection P > 0.01, < 3 beads in at least 5% of samples per probe, Non-CpG probes, single-nucleotide polymorphism (SNP)-related probes, multi-hit probes, and allosomal CpG probes were filtered out. In total, DNA methylation data were obtained for 763 541 probes. Background correction and beta mixture quantile normalization (BMIQ) to normalize the type I and type II probes was implemented using ChAMP (31). The singular value decomposition (SVD) method was used to assess batch effects and ComBat (32) was implemented for correction of multiple batch effects. Since genome-wide methylation analysis was performed on DNA obtained from buffy coats, their methylation status could potentially reflect the combination of blood cell types. RefbaseEWAS was implemented to correct for changes in distribution of white blood cells between different subpopulations using DNA methylation signatures in combination with a previously obtained external validation set consisting of signatures from purified leukocyte samples (28, 33). To reduce heteroscedasticity for highly methylated or unmethylated sites, β values were converted to M values in the lumi package (34) for further analysis calculated as (M = log2(β/(1−β))) (35). For association with maternal Hb levels, partial correlation with the aforementioned covariates was implemented.

Correlation of CpG methylation with gene expression

For correlation of CpG methylation with gene expression, data-frames were created for each gene tested with their respective CpG sites based on the annotation provided for the Illumina EPIC chip. The correlation of CpG methylation with gene expression was performed using Pearson correlation test in R after log transformation of expression CPM and methylation beta values. Gene-wise, false discovery rate (FDR) correction was applied by correcting the P values for the number of CpG sites tested for each gene.

Correlation of CpG methylation with phenotypes

Pearson correlation test was applied to assess the correlation of CpG methylation with phenotypes like placental weight, placental size, birth weight, waist-to-length ratio, cord insulin, and c-peptide levels after data normalization using log transformation. Only CpG sites whose methylation levels associated at least nominally with gene expression of their respective genes were selected for the analysis. Additionally, only genes showing at least nominal correlation with their respective trait were considered and FDR correction was applied by taking into account only the CpG sites tested.

Genome-Wide Association Study

DNA samples were genotyped on the Illumina Infinium OmniExpressExome-8 v1.4 using Illumina iScan at the Department of Clinical Sciences, Clinical Research Centre, Lund University, Malmö, Sweden.

Variants with call rate < 99%, extreme heterozygosity (> mean ± [3 × SD]), monomorphic variants or with low minor allele frequency (MAF < 1%) were removed before imputation. A cutoff of Hardy-Weinberg equilibrium (exact P < 5.7 × 10-8) and (exact was P < 10-4) was used for common and rare variants, respectively. Quality control was performed using the PLINK 1.9 software package (28). After quality control, 281 589 variants were left for imputation. The genome-wide association studies (GWAS) scaffold was mapped to NCBI build 37 of the human genome, and imputation to the 1000G reference panel (Phase 3- http://www.well.ox.ac.uk/~wrayner/tools/) for European ancestry SNPs were available after imputation. All SNPs with MAF < 0.01 and imputation quality < 0.4 were removed from analysis. Expression quantitative trait locus (eQTL) analysis was implemented in MATRIX-QTL after the stipulated transformation of data.

RNAseq and GWAS of Fetal Tissues

Human fetal pancreas (n = 16) and liver (n = 11) tissue was obtained from material available following elective termination of pregnancy during the first trimester at the University Hospital in Malmö, Sweden. The study was approved by the national ethics committee (permit number 2018-579), written and oral consent was obtained prior to collection. Tissue samples were dissected and homogenized in Trizol and further extracted using RNeasy columns for bulk extraction. The RNA integrity number (RIN) was > 7 and approximately 1 microgram of RNA was obtained upon Trizol extraction. RNA was purified and TruSeq Stranded Total RNA Ribo-Zero H/M/R Gold (Illumina) was used for library preparation. RNA sequencing was performed on the Hiseq 2000 (Illumina). Data were processed as per the GTEX pipeline (36). Briefly, the reads obtained were aligned to the human reference genome hg19/GRCh37 using STAR v2.4.2a based on GENCODE v19 annotation and gene level counts were obtained using RNA-SeQC v1.1.8 (37). The TPM values obtained were normalized and differential expression analysis was performed between embryonic and adult RNAseq data using edgeR, with focus on the specified genes whose expression was altered in UCB from early-anemic mothers.

Simultaneously GWAS was run on DNA extracted from the same embryos using the Illumina Omni Express Exome chip (Illumina). GWAS data was processed as described earlier and Matrix eQTL was run to identify genetic variation altering gene expression in each of the tissues (38).

Immunohistochemistry Assays

For immunochemistry, 6- to 8-month-old pancreata from human terminated fetuses were processed for paraffin (6-μm sections) and cryo-embedding (10-µm sections), respectively, and immunostained as described previously (39) with the following antibodies: guinea pig α-insulin (1:2000, Millipore/1:800, AB_433703), mouse α-glucagon (1:2000, AB_259852), mouse NUMBL (1:2000, AB_2851987), rabbit PIK3C2B (1:500, AB_2893430), rabbit P2RX7 (AB_2545496) antibodies. Nuclear counterstaining was performed using 4′,6-diamidino-2-phenylindole (DAPI; 1:6000, Invitrogen).

Gene Set Enrichment Analysis

Gene set enrichment analysis was performed using the Kolmogorov-Smirnoff tests to assess if the expression of a particular set of genes was enriched in EP anemia–exposed UCB compared with those from nonanemic mothers. For the gene sets, we tested for the enrichment of eGenes, or target genes from eQTLs in the UCB for SNPs previously associated with T2D, BW, and obesity. We also used the KEGG pathway gene lists to test for enrichment of hematopoietic pathways genes (40, 41).

Results

Expression Landscapes of Cord Blood From Offspring of Mothers With Anemia During Early Pregnancy Are Distinct From Controls

Mothers consistently having low Hb levels (Hb ≤ 10.9 g/dL) in early pregnancy (EP) defined as the first and second trimesters (21 weeks + 5 days) comprised EP anemic cases. These women were older, and had lower BMI (kg/m2) than control mothers without EP anemia (Fig. 1, Supplementary Table S1 (23)).

Figure 1.

Overview of study design and findings.

To examine differences in gene expression in UCB from newborns who had been exposed to these conditions, we performed RNA sequencing (RNAseq) from 48 UCB samples from newborns of mothers with who had EP anemia and 50 controls without EP anemia. We identified 137 differentially expressed genes, including 111 protein-coding genes and 26 non-coding and antisense RNAs including the LCORL, NMD3, and NMBP3 genes (Figs. 1, 2A, 2B, Supplementary Table S2 (23)). The majority of the genes were upregulated in UCB from newborns exposed to early maternal anemia.

Figure 2.

(A) Volcano plot showing differentially expressed genes (DEGs) between cord blood from mothers with early pregnancy anemia compared to controls. All genes showing FDR < 0.05 are presented in green. DEGs showing logFC > 1.5 are presented in orange. DEGs with FDR < 0.05 and log FC > 1 are presented in red while those with FDR < 0.005 and logFC > 1 are labeled. (B) The expression of LCORL, NMD3 and NCBP1 genes was upregulated in umbilical cord blood from offspring of mothers with early pregnancy anemia compared to controls. (C) cg12884187 and (D) cp06549407 DNAm negatively correlated with LCORL expression (E) cg12884187 DNAm correlated with maternal HB levels in early pregnancy (F) cg06549407 DNAm correlated positively with delivery length (G) LCORL expression correlated with delivery length.

We next examined if expression of these genes was also associated with maternal Hb levels at delivery. One gene, ZNF417, was significantly (FDR < 0.05) associated, whereas 29 other genes were nominally (P < 0.05) associated with maternal Hb at delivery, and the remaining 107 genes did not show any association (Supplementary Table S3 (23)).

In the search for pathways mediating the effect of anemia on gene expression in UCB, we performed pathway analysis using PANTHER(v14) (42). This analysis yielded an overrepresentation of genes involved in developmental, immune, and metabolic processes (Supplementary Figure S1a), including genes involved in integrin, cytokine and chemokine, endothelial, and P53 signaling, and apoptosis pathways (Supplementary Figure S1b (23)).

Epigenetic patterns accompany changes in gene expression in cord blood from mothers with EP anemia

EP anemia can modify fetal programming through epigenetic changes in the fetus, which in turn leads to changes in gene expression. We first aimed to identify CpG sites with DNA methylation (DNAm) associated with the expression of the 137 DEGs. The DNAm at 262 CpG sites in 41 genes were significantly associated with expression of their respective genes (FDR < 0.05). Some of these included cg11265882, cg12884187, and cg06549407 (Fig. 2C, 2D) which negatively correlated with LCORL expression, and CpG sites cg13580487 and cg08510057 in NMD3, which negatively correlated with NMD3 expression (Supplementary Table S4 (23)).

We then sought to investigate if the CpG DNAm which had correlated with gene expression here also associated with maternal Hb levels at the first contact (at GA 60.46 days [± 19.09]) in the first trimester. DNAm at these sites correlating with maternal Hb levels in the first trimester would indicate that epigenetic signatures that were altered due to early maternal anemia exposure and could be potential epigenetic signatures of fetal programming. When all CpG sites with nominal methylation-expression correlation were considered, DNAm at 42 CpG sites in 27 genes was also associated with maternal Hb levels in EP (Supplementary Table S5 (23)). Of particular interest was the CpG site cg12884187 whose methylation was negatively correlated with LCORL expression. Methylation at the same CpG site here correlated positively with maternal Hb levels in EP (Fig. 2E, Supplementary Table S5 (23)). DNAm at CpG sites in COBLL1, PIK3C2B, PRC1-AS1, P2RX7, FRY, and CIT was associated with gene expression levels of their respective genes and correlated with maternal Hb levels in EP (Table S5 (23)).

Genes Showing Altered Expression in EP Anemia–Exposed Cord Blood Correlates With Offspring Birth Weight

Maternal anemia in EP is associated with low BW (43, 44). Consistent with previous findings, maternal Hb levels at first contact (GA 60.46 [± 19.09] days) in our study were positively correlated with BW of the offspring (r = 0.26, P = 0.009) (Supplementary Figure S2A (23)) whereas maternal Hb at delivery was not (r = 0.17, P = 0.86) (Supplementary Figure S2B (23)).

We first investigated the correlation of LCORL methylation and expression with BW and placental weight and size, since: (a) LCORL showed differential expression (DE); (b) DNAm at CpG sites in LCORL associated with its expression; (c) DNAm at the same CpG site was also associated with maternal Hb levels in EP; and (d) LCORL is a known BW locus.

DNAm at 11 CpG sites correlated with LCORL expression (Supplementary Table S4 (23)). Of these, DNAm at 2 CpG sites, including cg10374588 and cg06549407 (Fig. 2F) (23), also correlated with length of the neonate. DNAm at CpG sites including cg11265882 and cg01535003 correlated significantly whereas DNAm at cg10374588, cg12884187, cg27176246, and cg01089331 correlated nominally with placental size (Supplementary Table S6) (23). DNAm at cg01535003 correlated significantly, whereas DNAm at cg06549407 and cg20358902 correlated nominally with placental weight (Supplementary Table S7 (23)). Expression of the LCORL gene was significantly associated with length of the neonate at birth (Fig. 2G), placental weight and size but not with BW (Supplementary Table S6 (23)).

We then tested the 137 DEGs for their correlation with BW; the expression of KLHDC1 significantly correlated with BW (Supplementary Figure S3A (23)) whereas FRY, MTA1, C2CD2L, RPS6KA5, and NCBP3 associated nominally (Supplementary Table S6 (23)). Methylation at one CpG site which previously correlated with MTA1 expression here nominally correlated with BW (Table S7 (23)).

Next, we then explored if any of these genes correlated with placental weight, placental size, and neonatal anthropometry. The expression of all 6 genes including KLHDC1 (Supplementary Figure S3B (23)) was also correlated with neonatal weight-length ratio, while that of KLHDC1, FRY, RPS6KA5, and NCBP3 correlated with abdominal circumference. FRY, C2CD2L, and NCBP3 correlated with head circumference while FRY and NCBP3 expression correlated with length of the neonate. In support of this, DNAm at cg07808960 and cg15212210 nominally correlated with weight to length ratio and head circumference, respectively.

Expression of 4 genes including EXOSC9, LCORL, PRC-AS1, and NMD3, which earlier showed differential expression in both studies, here correlated nominally with placental size as well as placental weight. FRY expression was also nominally correlated with placental size (Supplementary Table S6 (23)). Among the 9 CpG sites whose DNAm correlated with NMD3 expression, 3 CpG sites correlated with placental weight. Finally, CpG sites in C2CD2L and FBXL20 also correlated with placental size (Table S7 (23)).

Genes Showing Altered Expression in EP Anemia–Exposed Cord Blood Correlates With Insulin and C-Peptide Levels in Cord Blood

Insulin is key in regulating glucose homeostasis, and fetal insulin levels are also known to stimulate fetal growth (45). In our study, insulin levels were lower in UCB from EP anemic mothers compared with the controls, indicative of impaired insulin secretion among EP anemic offspring (Supplementary Table S1 (23)). We therefore sought to determine which of the genes previously consistently correlated with neonatal outcomes (BW, weight to length ratio, and others: KLHDC1, FRY, MTA1, C2CD2L, RPS6KA5, NCBP3, PRC1-AS1, and NMD3), also correlated with cord insulin levels.

Of the above listed 8 genes, expression of KLHDC1 significantly correlated with insulin levels (Supplementary Figure S3C) whereas NMD3 correlated nominally (Supplementary Table S8 (23)).

Assessing KLHDC1 expression in RNAseq data from 191 donor human pancreatic islets from another study, we found that KLHDC1 expression was not associated with that of insulin expression (Figure S3D) and was downregulated in T2D donor islets (Figure S3E).

The expression of NMD3 was nominally downregulated in T2D islets and negatively correlated with INS, the latter in both islets well as in cord blood (Table S8, S9) (23).

Genes Showing Altered Expression in EP Anemia–Exposed Cord Blood Associates With Insulin Secretion and Beta-Cell Function

Next, we examined if any of the 137 genes showed potential roles in insulin secretion and beta-cell function.

Look up of expression in RNAseq from 191 pancreatic donor islets

We first looked up the expression of all 137 DEGs in RNAseq data from islets from 191 donors (Islet GeneView). The expression of the COBLL1 gene was significantly downregulated in islets from T2D donors (n = 191) and correlated negatively with INS expression. SSTR5-AS1 and ELN expression was upregulated in T2D islets and correlated positively with INS (Fig. 3A-3F, Supplementary Table S9 (23)).

Figure 3.

RNAseq from human pancreatic islets (n = 188). COBLL1 (A), SSTR5-AS1 (B), and ELN (C) were differentially expressed between diabetic donor islets compared to controls. The expression of COBLL1 (D) was negatively whereas those of SSTR5-AS1 (E) and ELN (F) were positively correlated with INS expression. ZDHHC14 (G) expression was upregulated whereas LCORL (H), EZH1 (I) and DBF4B (J) expression was downregulated in T2D compared to non-T2D donor islets, ZDDHC14 (K) expression correlated positively with insulin expression whereas EZHI (M) and DBF4B (N) correlated negatively. LCORL (L) showed no correlation. Normoglycemic islets (n = 31) were exposed to normal (5.5 mmol/L) and high (18.9 mmol/L) for 24 hours. ZDHHC14(O) expression was significantly upregulated whereas that of LCORL (P), EZH1 (Q) and DBF4B (R) was downregulated upon high glucose stimulation.

The expression levels of 14 genes including COBLL1, BRCA1, FBXL20, PIK3C2B, SDHC, NMD3, EZH1, and DBF4B were negatively correlated whereas 5 genes including BRD9, NUMBL, and ZDHHC14 were positively correlated with INS expression (Fig. 3K-3N), Supplementary Table S9 (23)).

Look up of expression in RNAseq from glucose-treated pancreatic donor islets from hyperglycemic and normoglycemic donors

In order to test if any of these genes were also responsive to glucose, we then looked up these genes in our RNAseq data from normal (5.5 mmol/L) and high glucose (18.9 mmol/L) treated islets from normoglycemic donors and hyperglycemic donors from our previous study (46).

COBLL1 and SSTR5-AS1, whose expression was altered in T2D donor islets, here showed no response to acute glucose exposure in either hyperglycemic or normoglycemic donor islets. These could therefore be considered to have causal involvement in T2D pathogenesis.

Among the 14 genes that negatively correlated with INS expression, the expression of CATSPER2 and DBF4B was downregulated in 18.5 mM glucose-treated islets from normoglycemic donors. Among the genes that were positively correlated with INS, PLEKHJ1, and PRR26 were upregulated. Since these genes were expressed in response to acute glucose in the short term but not the long term, and they correlated with INS expression in islets, they could therefore be directly involved in insulin secretion (Supplementary Table S9 (23)).

Among the 12 genes showing nominal differential expression in T2D donor islets, 10 were also associated with INS expression. Among these, ZDHHC1 expression was upregulated whereas EZH1, LCORL, and DBF4B was downregulated upon glucose treatment (Fig. 4O-4R, Supplementary Table S9 (23)). These are potential key players in beta-cell function and likely to be important in the pathogenesis of T2D.

Figure 4.

(A) Expression of genes whose expression in altered in cord blood from mothers with early anemia in fetal vs adult pancreas. (B) Expression of selected genes in sorted fetal beta, alpha and adult beta-cells (C) Immunohistochemical staining of 8-week fetal pancreas was performed for P2XR7 (red), NUMBL (red), PIK3C2B (red), INS (green), and GCG (white). Scale bar indicates 50 µm, pictures were taken with a 20× objective. Arrows denote insulin positive cells. Note: False positive green staining outside pancreatic epithelium is due to auto-fluorescence from red blood cells.

Gene Expression in Embryonic Pancreas and Liver

Based on the hypothesis that genes whose expression was enriched in embryonic/ fetal tissue are more likely to be involved in development, we then examined the expression of these 137 DEGs in another study from Sweden where we had access to RNAseq data for fetal pancreas and liver samples. Here, we assessed their expression in embryonic pancreas and liver compared to expression in adult pancreas and liver respectively, using the publicly available GTEX database (36).

Pancreas

We next examined the expression of the 137 genes in 16 fetal (Swedish) and 328 adult pancreas (GTEx) to understand their role in development and/or function. Beta-cells in the pancreatic islets produce insulin. We reasoned that genes whose expression was high in fetal pancreas, and higher in fetal beta-cells compared to adult pancreas, and adult beta and fetal alpha cells respectively potentially reflect those that influence differentiation and ultimately beta-cell mass. Lower beta-cell mass is unable to compensate for the higher insulin demands due to insulin resistance, potentially due to obesity in later life, resulting in the development of diabetes (47, 48). To achieve this, we first performed differential analysis between fetal and adult pancreas to identify genes showing higher expression in the respective tissues. Next, we assessed their expression in human pancreatic islets, as well as sorted alpha and beta-cells from adult and fetal islets.

Discernable levels of expression in either the fetal or adult pancreas was determined for 121 out of the 137 genes. Among these, the expression of 43 genes was upregulated in the fetal relative to the adult pancreas, suggestive of their role in development, which included the LCORL, ATP5H, RPS6KA5, and other genes associated with neonatal outcomes, including FRY (Fig. 4A), Supplementary Table S10 (23)). Of these, the HELLPAR, ARHGAP20, and ZBTB8B genes showed fetal pancreas specific expression, showing ≥ 1 cpm in 75% of the samples in fetal pancreas but ≤ 1 cpm in at least 75% of the adult pancreas samples.

Expression in pancreatic donor islets

A total of 12 genes eg, ITGAD, KIAA2012, A2ML1, and P2RX7 showed very low to no expression in adult islets (43) (Supplementary Table S11 (23)). LCORL expression was downregulated in T2D islets (Fig. 3F).

Expression in sorted fetal and adult alpha and beta-cells

Next, we performed a selection of genes based on expression in RNAseq data fetal and adult sorted alpha and beta-cells reported by Blodgett et al (49) to focus on genes potentially involved in alpha and beta-cell development. We first compared our list of 43 genes upregulated in fetal pancreas with clusters that defined enrichment of fetal beta, alpha, or fetal in general. The , and TP53 genes were among a cluster of genes defined to be beta-cell specific whereas the LCORL was alpha-cell enriched. and MAP4K5 were fetal enriched (Fig. 4B).

Next, we assessed if our set of 43 genes showed higher expression fetal beta-cells compared to adult beta as well as adult alpha cells. The expression of 18 genes eg, A2ML1, FBXL20, NUMBL, P2RX7, PIK3C2B genes was higher in fetal beta compared to adult beta-cells (Blodgett et al (49), Table S12, Fig. 4B). Moreover, the expression of A2ML1, GPNMB, NUMBL, , PIK3C2B, ADAMTSL5, ITGAD, NLRP1 and PKD1L2 was higher in fetal beta compared to alpha cells (Fig. 4B, Supplementary Table S12 (23), Blodgett et al (49)).

Identification of genes potentially impacting beta-cell mass

In total, we found 15 genes based on: (1) higher expression in fetal pancreas compared to adult; (2) very low levels of expression in adult pancreatic islets; (3) being part of a cluster of genes defining fetal beta-cells; or (4) upregulated in fetal beta compared to adult beta-cells (Fig. 3B, Supplementary Table S13 (23)). Immunohistochemical staining of 8 weeks fetal pancreas showed that P2XR7 is mainly expressed in endocrine progenitors, NUMBL, in the whole pancreatic epithelium, including some insulin and glucagon producing cells, whereas PIK3C2B expression was detected in the epithelium (Fig. 3C).

Potential long-term impact on insulin secretion and T2D

We next examined if these genes had eQTLs in human pancreatic islets (n = 191) and found 3 genes to have genome-wide significant eQTLs (Supplementary Table S14a (23)). Next, we examined the loci spanning the genes for genetic associations with indices of insulin secretion and T2D in the AMP knowledge portal. Of these, the rs34815312 SNP, which is an eQTL for the P2RX7 gene in islets, was associated with insulin secretion indices in a previous intravenous glucose tolerance test GWAS study (P = 0.03), with T2D in FinnGen GWAS (P = 0.01), HbA1c in TOPMed meta-analysis (P = 0.01) and a significant association with red blood cell count in the PheWAS scan of the UK biobank (P = 0.00002). The eQTLs for A2ML1 associated with type 1 diabetes and T2D from previous GWAS (Supplementary Table S13 (23)). Multiple nominal eQTLS for these genes also associated with insulin secretion indices in the MAGIC dataset (Supplementary Table S14b (23)).

We suggest that these genes can potentially lower beta-cell mass, which may underlie lower insulin secretion in response to higher demands for insulin in later life, thus leading to development of T2D.

Liver

Fetal liver is the site of hematopoiesis in mammals. Therefore, we compared RNA expression data from 11 human embryonic/fetal (7-14 weeks postconception) and 226 adult liver (GTEx database), and observed that for the 137 above genes, 39 genes, for example, TP53, RPS6KA5, and PRC-AS1, were expressed in fetal liver, but showed low to no expression in adult liver. In contrast, 23 genes showed higher expression in adult liver compared to fetal liver (Supplementary Table S15, Figure S4A (23)). We then compared the expression of the 39 genes showing higher expression in fetal liver to gene profiles that modulated specification of erythroid, megakaryocyte and mast cell lineages (50). TARBP1, TP53, RPS29, ATP5E, and EXOSC9 showed almost ubiquitous expression, whereas ATP6V0A2 expression was seen in mast cells and megakaryocytes, CHST12 in megakaryocytes natural killer (NK) and endothelial cells, RPS6KA5 in mast cells, CBX5 in early erythrocyte, Hematopoietic Stem Cell-Multi-Potent Progenitor cells (HSC_MPP), mid-erythroid, NK, pre-B, and megakaryocytes, and RAD50 in mid-erythroid cells (Supplementary Figure S4B, S4C (23)).

T2D-Associated Genes Overrepresented Among Genes Expressed in Cord Blood From Mothers With Early Anemia

Poor in utero conditions have been associated with long-term outcomes in the offspring for eg, T2D, obesity, and retinopathy. In order to evaluate if maternal EP anemia alters expression of genes associated with long-term outcomes, we first identified if the known T2D, obesity, and hematopoietic genes had eQTLs in the UCB and if their expression was impacted due to maternal early anemia.

Previous GWAS have identified 403 loci associated with T2D risk (Supplementary Table S16 (23)). Among them, 252 SNPs showed nominal cis-eQTL effects on 359 genes in UCB (Supplementary Table S17 (23)). Of them, 72 genes showed high expression in fetal compared to adult pancreas (Supplementary Table S18 (23)). Gene set enrichment analysis showed that these 72 genes were significantly overrepresented in UCB from mothers with EP anemia (P < 0.021, Supplementary Figure S5 (23)).

Gene set enrichment analysis using eGenes in UCB from previous GWAS associated BW and obesity loci as well as genes involved in hematopoiesis from KEGG pathways did not show any significant enrichment (Supplementary Figures S6-S9 (23)).

Discussion

By leveraging global gene expression, DNA methylation, and genotyping on UCB obtained from offspring of mothers who had anemia during EP, we provide novel insights into fetal programming signatures and molecular pathways related to metabolism and T2D. Maternal anemia during EP is associated with low BW, small for gestational age, and preterm delivery (12, 51, 52). Through differential expression analysis, we found the expression of 137 genes altered and 8 of these changes were also observed in our second study. The expression of some of these genes also associated with maternal Hb levels at delivery, and it is possible that these expression changes initiated in EP are retained throughout pregnancy and are also modulated by maternal Hb levels at a later stage.

The developmental plasticity in response to nutritional depletion in EP due to maternal anemia is speculated to alter fetal programming and cause adverse consequences in later life (53). This was suggested to be mediated by epigenetic changes in developmental genes (53). We here also show accompanying epigenetic signatures that correlated with gene expression that are potentially altered in UCB of offspring due to early maternal anemia in both the cohorts.

In accordance with previous reports showing an association between anemia and small for gestational age and low BW (43, 51, 54, 55), we observed an association between maternal Hb levels in early anemia and BW of the offspring. We here provide data showing how epigenetic mechanisms could be involved in this scenario. The Ligand Dependent Nuclear Receptor Corepressor Like (LCORL) gene is a component of the polycomb repression complex-2 (PRC2), which mediates methylation of histone H3lys27 (H3K27), and is important for determining methylation status and regulation of cellular identity during fetal development (56, 57). DNA methylation at CpG sites in the LCORL gene was robustly associated with LCORL expression, and also with maternal Hb levels in EP. The epigenetic signature is a potential consequence of maternal anemia in EP and could therefore represent an example of a signature of fetal programming. Moreover, LCORL expression was associated with placental size and weight, which in turn could affect fetal growth and BW (58-60). Furthermore, CpG methylation correlated with LCORL expression, and also correlated with length of the newborn, placental size, and placental weight. This provides further evidence for the mechanism by which altered fetal programming could affect the outcomes for the newborn with potential for long-term consequences.

Interestingly, insulin levels were lower in UCB from neonates of EP anemic mothers compared to those without. An explanation for this could be that undernutrition associated with EP anemia evokes an impaired fetal pancreatic development and maturation resulting in reduced in vivo insulin secretion and subsequently increased lifetime risk of T2D. Insulin is furthermore an important growth hormone, which has a significant effect on fetal growth, mediated through altered gene expression. KLHDC1 represents such an example, whose expression was increased in UCB from early-anemic mothers, and was negatively correlated with BW, weight to length ratio of the neonate, and with UCB insulin levels. In adult pancreatic islets, KLHDC1 expression correlated negatively with BMI and did not associate with insulin expression. We hypothesized that higher levels of KLHDC1 expression could be a result of lower insulin and therefore result in poor fetal growth leading to lower BW and weight to length ratio. While the KLHDC1 gene is expressed in embryonic and adult liver and pancreas, its expression is markedly higher in the adult tissues. The KLHDC1 gene encodes the Kelch domain containing 1 protein, which is highly expressed in pancreas, liver, and skeletal muscle where it is suggested to have roles in skeletal muscle development (61). Low BW has been clearly associated with higher risk of T2D and CVD (62). We hypothesized that KLHDC1 is involved in fetal programming induced by early anemia during pregnancy that alters future risk of metabolic disease through BW and fetal development, and this is potentially mediated by insulin.

Indeed, we found many of the genes whose expressions were altered in UCB from newborns of EP anemic mothers to also correlate with insulin levels in human pancreatic islets, and some whose expression was also altered in islets from T2D donors. Given the importance of insulin as a growth factor, we can therefore speculate that these genes affect fetal growth through altered insulin secretion.

Many of the DEGs have established roles in development. In support of this, multiple genes whose expression was altered in UCB with anemia exposure during EP also showed elevated expression in human embryonic pancreas and liver and low expression in their adult tissues. Of particular interest were genes potentially involved in regulating beta-cell mass. We show how some genes showing altered expression in UCB of offspring with early maternal anemia are also expressed in fetal pancreas, but show low to no expression in adults, and are enriched in fetal beta-cells, and associating with fetal insulin levels. These are the most likely candidates to influence beta-cell mass. People with T2D have deficits of beta-cell mass ranging from 20% to 65% (63, 64). While beta-cell death might explain this deficit, it has been hypothesized that insufficient islet development during the prenatal and postnatal period could play a key role. While this might not manifest during early adulthood, insufficient insulin secretion in later life in response to developing obesity and increased insulin demands lead to the development of T2D (65). Indeed, some of these dysregulated loci including the P2RX7, KIAA2012, A2ML1, and other genes potentially involved in beta-cell development also associated with indices of insulin secretion and T2D risk in GWAS studies (66, 67). Moreover, some genes have previously been shown to have developmental roles. NUMB Like Endocytic Adaptor Protein coding NUMBL is an important gene for cell differentiation and embryonic neurogenesis by inhibiting NOTCH signaling (68). ZDHHC14, a protein-cysteine S-palmitoyltransferase coding gene, is a known tumor suppressor. PIK3C2B plays significant roles in signaling pathways involved in cell proliferation, oncogenic transformation, cell survival, cell migration, and intracellular protein trafficking (69). This lower beta-cell mass could potentially explain the lower cord insulin levels in offspring of EP anemic mothers.

The liver is the site of hematopoiesis in the mammalian embryo up until the middle of the second trimester (50) and some of our putative fetal programming genes were also expressed in precursors of hematopoiesis, although none of them were key players reported in Popescu et al (50).

Hb is the key molecular transporter of oxygen from lungs to all tissues. Under normal conditions, the amount of oxygen delivered far exceeds the demand. During pregnancy, oxygen demands of the fetus are substantially high and hypoxia-inducible factor-1α (HIF-1α) protein has been shown to accumulate in placentae and hearts of anemic fetuses of animal models (70, 71). The ATP5H gene coding the subunit D of mitochondrial H(+)-ATP synthase, a multi-subunit complex, and ADCY9 gene encoding the enzyme adenylate cyclase type 9 are genes important in ATP formation, and both are upregulated in early maternal anemia–exposed UCB. These genes are vital for mitochondrial energy production (72). Therefore, we speculated that ATP5H upregulation is a consequence of hypoxia subsequent to anemia during EP potentially as a compensatory mechanism to combat hypoxia. In addition, cyclic AMP (cAMP) is an important player in the regulatory mechanism of hypoxia response by downregulating the activity of the master regulator of hypoxic activation of genes for angiogenesis, hormone synthesis, glycolysis, and cell survival, the hypoxia-inducible factor 1 (HIF1) gene (73). Interestingly, ADCY9 belongs to the protein kinase A signaling pathway and was deregulated leading to transgenerational impairment of ovarian development in fish (74).

In the absence of data on development of metabolic disease in the offspring, we examined expression levels in the UCB samples of validated T2D susceptibility genes (derived from eQTL SNPs of T2D loci) (66). Indeed, we observed significant enrichment of T2D-associated genes among the genes that in UCB samples were differentially expressed in cord blood from neonates exposed to EP maternal anemia. This provides compelling evidence of a potential overlap between mechanisms involved in developmental programming of T2D by anemia in pregnancy, on one side, and genetic mechanisms underlying T2D susceptibility on the other. In further support of this idea, the expression of some of these T2D susceptibility genes was altered in pancreatic islets from T2D donors.

Our studies were focused on –omics signatures in UCB, since UCB is representative of the intrauterine environment and relatively easily accessible. Moreover, UCB DNAm and changes in gene expression have been associated with several outcomes linked to early childhood weight and adiposity (75) and autism spectrum disorder (76). Finally, only few whole transcriptomic studies on UCB have been performed to date, and none of these have taken maternal EP anemia into account (77).

This study has potential limitations. One such limitation is the limited power and lack of precise replication cohorts. Another limitation is the lack of long-term outcome data in the offspring. This limitation, however, we share with most other studies, particularly those based in developing countries in the field. As a surrogate, we use BW as a predictor of long-term outcome as previously robustly established. Altered expression of important developmental genes, like imprinted genes, could have potential effects on long-term outcomes. Interestingly, the expression of T2D-associated genes was enriched in UCB from anemic mothers. A potential limitation is that most of our data are from UCB and it is not entirely clear how well UCB reflects processes in other organs.

We cannot in the current study design rule out any association between anemia diagnosis in the pregnant women, and whether the women also were subjected to undernutrition. Indeed, the association between undernutrition and fetal programming mechanisms is one of the most studied within the Developmental Origins of Health and Disease hypothesis. Anemia is highly linked to low iron intake, and that can possibly be associated with malnutrition. Several studies have investigated the relationship between maternal undernutrition and low BW, differential expression, and methylation in offspring. However, in our thorough literature search we did not detect any publications that described association between maternal undernutrition and the identified gene candidates (LCORL, P2RX7, PIK3C2B, and NUMBL).

The women participating in the study were all given offered daily iron and folic acid supplementation from pregnancy recognition till delivery. If anemia was diagnosed, it was treated with increased supplementation with either 2 to 3 combination tablets of iron-folic acid per day or with Hemovit multivitamin syrup (200 mg ferrous sulfate, 0.5 mg vitamin B6, 50 µg B12, 1.5 mg folic acid, and 2.33 mg zinc sulfate per 5 mL [Shelys Pharmaceuticals, Dar es Salaam, Tanzania]) in a dose of 10 mL 2 to 3 times daily depending on severity (22). Therefore, this could be a potential confounder and its effects in gene expression and DNA methylation cannot be ruled out.

In conclusion, we provide a catalog of effects on how maternal early anemia in pregnancy causes changes in the expression of key genes and affects BW with potential impact on future metabolic disease. The study was carried out in Korogwe, Tanzania, where a high incidence of anemia and infectious diseases in the parent generation is followed by higher rates of noncommunicable diseases in subsequent generations (78). The intrauterine period represents a key trigger of later events and therefore provides a window of opportunity for interventions, aiming at reducing untoward effects in the offspring.