Fetal HLA-G mediated immune tolerance and interferon response in preeclampsia.

BACKGROUND: Fetal immune tolerance is crucial for pregnancy success. We studied the link between preeclampsia, a severe pregnancy disorder with uncertain pathogenesis, and fetal human leukocyte antigen G (HLA-G) and other genes regulating maternal immune responses.
METHODS: We assessed sex ratios and regulatory HLA-G haplotypes in population cohorts and series of preeclampsia and stillbirth. We studied placental mRNA expression of 136 genes by sequencing and HLA-G and interferon alpha (IFNα) protein expression by immunohistochemistry.
FINDINGS: We found underrepresentation of males in preeclamptic births, especially those delivered preterm or small for gestational age. Balancing selection at HLA-G associated with the sex ratio, stillbirth, and preeclampsia. We observed downregulation of HLA-G, its receptors, and many other tolerogenic genes, and marked upregulation of IFNA1 in preeclamptic placentas.
INTERPRETATION: These findings indicate that an evolutionary trade-off between immune tolerance and protection against infections at the maternal-fetal interface promotes genetic diversity in fetal HLA-G, thereby affecting survival, preeclampsia, and sex ratio. We highlight IFNA1 as a potential mediator of preeclampsia and a target for therapeutic trials. FUNDING: Finnish Medical Foundation, Päivikki and Sakari Sohlberg Foundation, Karolinska Institutet Research Foundation, Scandinavia-Japan Sasakawa Foundation, Japan Eye Bank Association, Astellas Foundation for Research on Metabolic Disorders, Japan Society for the Promotion of Science, Knut and Alice Wallenberg Foundation, Swedish Research Council, Medical Society Liv och Hälsa, Sigrid Jusélius Foundation, Helsinki University Hospital and University of Helsinki, Jane and Aatos Erkko Foundation, Academy of Finland, Finska Läkaresällskapet, Novo Nordisk Foundation, Finnish Foundation for Pediatric Research, and Emil Aaltonen Foundation.

Evidence before this study

Preeclampsia, a hypertensive pregnancy disorder, is a major cause of maternal and perinatal mortality and morbidity. Its root placental etiology remains unclear. Tens of studies, with still inconclusive results, have searched for the link between preeclampsia and HLA-G, a shield antigen to prevent killing of fetal trophoblasts by maternal immune cells. We hypothesized that the proposed but still disputed sex bias in preeclamptic births might reflect failure in fetal HLA-G mediated immune tolerance in preeclampsia, leading to stronger maternal immune responses towards male fetuses, and demonstrate the particular role of HLA-G in human pregnancy.

Added value of this study

We studied a large population and clinical cohorts and demonstrated a deficiency of male births in preeclampsia. We linked fetal regulatory HLA-G haplotypes to the offspring sex and pregnancy complications through the major mechanisms of balancing selection: negative frequency-dependent selection and heterozygote advantage. We found sex-linked downregulation of HLA-G, downregulation of HLA-G receptors and many other tolerogenic genes, and increased interferon signaling in preeclamptic placentas. Our results suggest that fetal HLA-G modulates fetal survival, preeclampsia, and human birth sex ratio.

Implications of all the available evidence

Although Fisher's principle, a theoretical argument central in evolutionary biology, posits that negative frequency-dependent selection keeps sex ratios balanced, there is very limited previous evidence across all species for this selection affecting any loci. This study suggests that balancing selection at the fetal HLA-G locus is associated with the human birth sex ratio and with a common, poorly understood and managed pregnancy disorder, preeclampsia. Our results provide support for immune-mediated mechanisms of preeclampsia and suggest that interferon signaling could be targeted in treatment and prevention of the disease as shown in special cases of preeclampsia.

Alt-text: Unlabelled box

Introduction

Whether human male-to-female sex ratios exhibit appreciable variation has been a subject of debate in the biological sciences. Recent large-scale data have confirmed an unbiased sex ratio at conception: equal numbers of X-bearing or Y-bearing sperm fertilize human eggs. However, the global birth sex ratio of 106 males for every 100 females, already seen in 20-week fetuses, indicates that more females are lost in early human pregnancies. In contrast, male fetuses might be at higher risk for late miscarriages and stillbirths [1]. Together, these findings suggest that undetermined fetal sex-specific mechanisms contribute to human pregnancy success and phenotypic variation in the sex ratio arises during pregnancy.

In humans, balancing selection at the major histocompatibility complex (MHC) locus has been proposed to modulate a deficiency of human leukocyte antigen (HLA) homozygotes through maternal-fetal interaction, [2] without regard for sex-specificity. Increased HLA similarity in couples is associated with recurrent miscarriages, further supporting the need for the fetus to differ from maternal HLA [3]. Although the exact modifier locus remains unclear, only a limited pattern of HLA antigens (HLA-C, HLA-E, HLA-F and HLA-G) are expressed by fetal trophoblasts to prevent rejection by maternal immune cells [4]. Of these, the non-polymorphic HLA-G is the most studied due to its trophoblast-restricted expression and multiple isoforms inducing both local and systemic immunotolerance during pregnancy [5].

Tens of studies, with still inconclusive results, have searched for the link between HLA-G, miscarriages, and the hypertensive pregnancy disorder preeclampsia. While reduced HLA-G expression observed in the preeclamptic placenta is expected to facilitate maternal immune reactions to fetal alloantigens, its role in the pathogenesis is not yet clear [6], [7], [8], [9], [10], [11]. Maternal immune maladaptation in preeclampsia, as one of the major theories, is supported by the shift from antibody-mediated T-helper 2 (Th2) and regulatory T cells (Treg) to cell-mediated T-helper 1 (Th1) responses, and by aberrant natural killer (NK) cell activity [12]. While current pathogenetic models emphasize the heterogeneity of preeclampsia, rather than immunological mechanisms, [12,13] a possible sex bias in preeclampsia fetuses should turn attention toward fetal antigens. Remarkably, preeclampsia placentas from male fetuses show more inflammation, hypoxia, and apoptosis, [14] even if a recent meta-analysis involving over 3.1 million women failed to show association between male fetal sex and preeclampsia [15]. In contrast, several cohort and epidemiological studies have demonstrated an unexplained female predominance in early preeclampsia [16], [17], [18].

This study was motivated by the proposed but still disputed sex bias in preeclamptic births, and particularly, the low male/female ratio in children of both women and men descending from a proband with preeclampsia [19]. We hypothesized that male fetal loss might reflect failure in fetal HLA-G mediated immunotolerance in preeclampsia, stronger maternal immune responses towards male fetuses, and demonstrate the particular role of HLA-G in maternal-fetal interaction that affects fetal survival in human pregnancy.

Materials and methods

Birth registry data. The Finnish Medical Birth Registry, hold by the Finnish Institute for Health and Welfare in Finland, involves all births after 22 weeks of gestation and data on maternal and fetal characteristics. Data spanning 30 years, 1987 to 2016, and including information on all live births and stillbirths from 22 gestational weeks or from ≥500 g of weight, data on gestational ages at birth, and fetal sex were collected from both singletons and multiple pregnancies (n = 1,790,123). Registry-based diagnoses of mild, severe and unspecified preeclampsia, as recorded during the hospitalization for labor and delivery, were used to build up the birth cohort with preeclampsia (n = 38,752) and the remaining control population (n = 1,751,371) without preeclampsia. Registry-based data on previous maternal miscarriages were available from 36,340 (94%) of preeclamptic births and were assessed against 1,755,203 births without preeclampsia. In this dataset, the prevalence of miscarriages before preeclamptic births was 17.8% (n = 6,469), and 20.3% in the control population.

1000 Genomes data and HLA-G 3′UTR haplotype analysis. We used Ensembl [20] on the reference genome GRCh38 to access 1000 Genomes Project data (1000 Genomes: A Deep Catalog of Human Genetic Variation, RRID:SCR_006828) [21]. Genotypes for all HLA-G 3′ untranslated region (UTR) SNPs were retrieved and automated scripts were created to build up HLA-G 3′UTR haplotypes presenting global frequencies of >1% (UTR-1 to UTR-7, UTR-10, UTR-18; Fig. S1) [22]. Haplotypes were analysed and the sex ratios (male/female) of individuals carrying each combination of two haplotypes, diplotypes, were calculated. Of the total 2,504 individuals (male/female ratio 0.97), only ten individuals (0.4%) carried one or two rare HLA-G 3′UTR haplotypes and were not included in the analyses.

Stillbirth cohort and controls. Samples from the Swedish Biobank Karolinska involved DNA from stillborn fetuses (n = 277) [23] and unselected population controls (n = 622). The sex of the individuals was confirmed by PCR. We used validated primer pairs for the Y chromosomal gene SRY and for the X chromosomal gene ATL1 [24,25]. PCR assays were run in a volume of 15 μL containing 1 x PCR buffer, 200 μM of dNTPs, 0.1 μM of SRY primers, 0.15 μM of ATL1 primers, 1.0 U of HotStar Taq Plus DNA Polymerase (Qiagen, Valencia, CA), and 50 ng of genomic DNA. Cycling conditions were initial 95 °C for 5 min to activate the polymerase, followed by 35 cycles of 95 °C for 30 sec, 57.1 °C for 30 sec, and 72 °C for 60 sec, and by the final extension of 72 °C for 5 min. The PCR products were electrophoresed on 2% agarose gels. For the male-biased sampling of Swedish controls, we did not study the frequencies of diplotypes and the sex ratios between the stillborn fetuses (sex ratio 1.18) and the controls (sex ratio 1.21).

FINNPEC series and samples. The Finnish Genetics of Preeclampsia Consortium (FINNPEC) cohort involves over 1450 women having a singleton pregnancy with preeclampsia, defined as hypertension and proteinuria occurring after 20 weeks of gestation. Hypertension was defined as systolic blood pressure ≥140 mmHg or diastolic ≥90 mmHg. Proteinuria was defined as the urinary excretion of ≥0.3 g protein in a 24 hr specimen, or 0.3 g/L, or two ≥1+ readings on a dipstick in a random urine determination without evidence of a urinary tract infection. After the recruitment of a preeclamptic woman, a woman with a non-preeclamptic pregnancy attending the same clinic was recruited as a control. Women with chronic hypertension and superimposed preeclampsia, fulfilling the above-mentioned criteria of hypertension and proteinuria, were included among preeclamptic cases. Characteristics of the cohort, classification of mild and severe preeclampsia, and details on biological samples have been described [26]. In the FINNPEC cohort, there is only one pregnancy record for each woman, without repeat pregnancies. We utilized clinical data on previous pregnancies and offspring sex ratio (male/female). The numbers of the males and females at birth were tested based on the number of previous pregnancies (zero to three) among 1,538 offspring from preeclamptic births, and separately in 333 small for gestational age (SGA) offspring, in comparison with the 1.79 million births from the Finnish Birth Registry data. The sex ratio of the whole FINNPEC series of newborns (n = 1,247) was 1.04, while preeclamptic newborns showed the ratio of 0.92 and controls the ratio of 1.16.

DNA samples from newborns of preeclamptic (n = 586) and control (n = 661) pregnancies were analysed. For RNA extraction, we used 182 deep frozen placental biopsy samples, of which 163 were included in the final analyses. The biopsies were dissected samples of villous tissue under the basal plate [26]. Formalin-fixed and paraffin-embedded sections from 24 placentas (13 preeclampsia and 11 controls), one sample (n = 3) or two samples (n = 21) from each placenta taken shortly after cesarian section, were studied by immunohistochemistry. Third trimester maternal serum samples (n = 209; 158 preeclampsia and 51 controls), taken at variable time points during hospitalization (0 to 31 days before the delivery), and fetal cord plasma samples (n = 176; 66 preeclampsia and 110 controls), taken immediately after delivery, were used for ELISA.

Sequencing and HLA-G 3′UTR haplotype analysis. Haplotypes of HLA-G 3′UTR were determined by Sanger sequencing among newborns of the FINNPEC series and in the stillbirth cohort and controls using the previously described primers [27]. PCR assays were run in a volume of 10 μL involving 1 x Phusion HF buffer, 200 μM of each dNTPs, 0.5 μM of each primer, 0.02 U/μL of Phusion DNA-polymerase (Thermo Fisher Scientific, Waltham, MA), and 20 ng of genomic DNA. Cycling conditions were 98 °C for 30 sec, followed by 35 cycles of 98 °C for 10 sec, 64 °C for 30 sec, and 72 °C for 15 sec, and final extension of 72 °C for 7 min. ExoSAP-IT PCR clean-up was performed according to manufacturer's instructions (Thermo Fisher Scientific). PCR products were sequenced with the reverse primer, for the detection of the 14 bp deletion/insertion, utilizing the BigDye Terminator v3.1 Cycle Sequencing Kit and the Applied Biosystems 3730 Genetic Analyzer (RRID:SCR_018052). HLA-G 3′UTR haplotypes were assessed manually by two readers independently. Haplotypes presenting global frequencies of >1% (UTR-1 to UTR-7, UTR-10, UTR-18) [22] and the sex ratios (male/female) of individuals carrying each diplotype were assessed.

RNA extraction. Placental samples, stored snap frozen or in RNAlater (Thermo Fisher Scientific) in −80 °C, were thawed and homogenized using the FastPrep-24 5 G Classic Instrument (2 × 30 s at 4.5 m/s) and Lysing Matrix D tubes according to manufacturer's instructions (MP Biomedicals, Santa Ana, CA). Samples were lysed using 600 μl RTL buffer and 6 μl ß-mercaptoethanol on 20 mg of tissue homogenate, followed by RNA extraction using the RNeasy Mini Kit according to manufacturer's instructions (Qiagen). RNA concentrations were measured (Thermo Scientific Nanodrop 2000 microvolume spectrophotometer, RRID:SCR_018042) and samples diluted in 10 ng/μl. The 2100 Bioanalyzer Instrument (RRID:SCR_018043) (Agilent Technologies, Santa Clara, CA) was utilized for measurement of RNA integrity (RIN) values. For most samples, the RIN values ranged from 7 to 8.8, but samples with the lower range of 4 to 6.9 gave highly comparable results in targeted allele counting by sequencing (TAC-seq) assay testing, and were included in the final assays. RNA samples were quantified (Thermo Fisher Qubit 2.0 Fluorometer, RRID:SCR_018095) using Invitrogen RNA HS Assay Kit (5 to 100 ng) or Invitrogen RNA BR Assay Kit (20 to 1000 ng) according to manufacturer's instructions (Thermo Fisher Scientific).

Selection of genes and probe design for TAC-seq. In addition to the main target HLA-G, we selected several genes known to regulate HLA-G expression, [5] trophoblast markers, [28], [29], [30] markers of leukocytes selected based on an in-house data and FANTOM data (FANTOM DB, RRID:SCR_002678), [31] members of the B7 family, [32] imprinted placental genes, [33] preeclampsia-associated genes, [34], [35], [36], [37], [38], [39] and housekeeping genes (GAPDH, UBC, TOP1, TBP, SDHA, MRPS18A, and YWHAZ), in total number of 146 targets. After adding 14 additional targets to detect seven different isoforms of HLA-G [5], altogether 160 targets were included in probe design (Table S1).

TAC-seq probes were designed to target a 54 bp of coding sequence near the 5′ end of the genes using the software produced by the researchers of this study (https://hindrek.shinyapps.io/probe_design/). GC content of the probes was set preferentially at 45–55%. For HLA-G isoforms, the probes were designed manually to cover the splice sites characteristic of each isoforms and were tested not to overlap with other HLA I genes by the BLAST program [40].

TAC-seq library preparation. The selected genes were analysed using the TAC-seq method [41]. Each target including HLA-G (n = 146), alternative targets (n = 14) to detect different HLA-G isoforms, and ERCC spike-in probes (n = 7; ERCC-00004; 00136; 00108; 00116; 00092; 00095; 00131; Thermo Fisher Scientific) had a pair of specific oligonucleotides (Table S1). Total-RNA of 2 µl (50 ng/µl) was mixed with 2 µl denaturation buffer and heated 1 min at 80°C. After cDNA synthesis in 5 µl, 1 µl of 10 µM TAC-seq probe mixture was added and incubated 1 hr at 60 °C before adding thermostable ligase. PCR was carried out in a volume of 40 µl containing 1 × proofreading HOT FIREPol Blend Master Mix (Solis BioDyne, Tartu, Estonia) and 250 nM primers. Two-step size-selection was performed to get rid of abundant cDNA (>500 bp) and 80 bp TAC-seq specific by-product. The 180 bp library was quantified and visualized on TapeStation High Sensitivity D1000 ScreenTape (Agilent Technologies). TAC-seq libraries were sequenced with Illumina NextSeq 500 (RRID:SCR_014983) high output 75 cycles kit (Illumina, San Diego, CA) using 62 bp read-1 and 6 bp index reads.

TAC-seq data processing. The data were deposited in NCBI's Gene Expression Omnibus (GEO: GSE125460; RRID:SCR_005012) [42]. Raw data (BCL files) were converted to FASTQ files using the bcl2fastq2 Conversion Software (version 2.20; RRID:SCR_015058), allowing no mismatches in the index reads. The sequencing quality was evaluated using MultiQC (version 1.3) [43]. Quality controlled samples were analysed using the TAC-seq data analysis software (available online at https://github.com/cchtEE/TAC-seq-data-analysis) as described [41] allowing up to three mismatches per targeted region of a sequence. For absolute molecule counts, each read incorporated 4 + 4-bp unique molecular identifier (UMI) sequence [44]. To reduce the potential over-estimation of molecules due to the accumulated sequencing errors, only reads that appeared at least twice were counted as a unique molecule.

Starting from a total of 239 samples including replicates, four samples were removed as outliers based on their ERCC spike-in molecule counts, and ten samples were excluded due to their low total molecule counts of <300,000. Before the analyses, 57 replicate samples were excluded, and four individual samples showing low correlation with other samples were removed as outliers. Finally, one sample without successful HLA-G 3′UTR haplotype analysis was excluded. As a result, 163 individual samples were retained for further analyses: 81 severe preeclampsia, 19 mild preeclampsia, and 63 control samples. Clinical characteristics related to the samples and subjects are presented in Table S2.

As described above, we started with 167 targets, including seven spike-in probes, for TAC-seq analysis (Table S1). First, two targets (IL2, MMP3) were excluded due to low expression (no count in more than half of the samples). Next, 8 targets (B2M, GAPDH, HSD3B1, KISS1, PAGE4, PLAC4, SEMA3B, SERPINE2) were removed since they were considered to reach saturation. We measured HLA-G expression by 15 different probe pairs and found highly similar patterns of differential expression; also for the major soluble isoform (Fig. S2). To avoid bias for many probes for HLA-G, only one HLA-G target was included in the final analysis. The probes for this target (target name HLA-G; Table S1) were designed by the software for the unique sequence in exon 2, which is present in all isoforms, and those 14 other targets for different HLA-G isoforms were excluded. Finally, probes for 136 targets retained for the data analysis.

We compared HLA-G mRNA expression levels, adjusted by offspring sex, mode of delivery, and gestational age at birth, in preeclampsia and controls by using the Welch's t-test, and accordingly, studied expression levels related to different HLA-G 3′UTR diplotypes using linear regression. Moreover, comparison of placental HLA-G expression with 3′UTR haplotypes was performed in two subgroups, based on the linear regression analysis and residuals of the observed and predicted sex ratios. The haplotype combinations with positive residuals were defined as male excess and those with negative residuals as female excess. The combination of UTR-3/UTR-5, for which three males and no females were found, was defined as male excess but was not shown in the graph. The mean expression levels of HLA-G mRNA from the unadjusted TAC-seq sequencing output were compared between the two haplotype groups, male and female excess, by the Welch's t-test in preeclamptic (mild and severe) and control placentas.

Differential expression analysis between control and severe preeclampsia samples was performed with the ‘nbinomWaldTest’ method of an R (version 3.5.0) package DESeq2 (version 1.20.0; RRID:SCR_015687) [45] using the following housekeeping genes as ‘controlGenes’ for ‘estimateSizeFactors’: UBC, TOP1, TBP, SDHA, MRPS18A, and YWHAZ. Here, the offspring sex, mode of the delivery (cesarian section or vaginal birth), gestational age at birth (divided in three stages: <32 weeks, 32–36 weeks, ≥37 weeks), and disease status (control, mild preeclampsia, severe preeclampsia) were included in the design model to correct the potential confounding factors. IGFBP1, ITGB5, KRT6A, MMP1, MMP9, and TGFB1 were regarded as outliers due to the Cook's distance >0.3 and excluded from the differential expression analysis. A gene was considered as differentially expressed when the Benjamini–Hochberg adjusted P-value was <0.05.

The expression data were normalized with DESeq2 (RRID:SCR_015687) and transformed into log2 (x + 1). These values were then adjusted by the offspring sex, mode of delivery, and gestational age at birth as described above using the ‘removeBatchEffect’ function of the LIMMA (version 3.36.5; RRID:SCR_010943) [46].

The plots were generated with ggplot2 (version 3.0.0; RRID:SCR_014601). The ROC curve (AUC) was calculated using the pROC package (version 1.13.0) [47]. The statistical significances of the differences between two ROC curves were tested by bootstrap method with 2,000 replicates using the ’roc.test’ function.

STRT RNA-seq library preparation and data processing. We validated the TAC-seq experiment by STRT RNA-seq by studying RNA samples from 32 preeclamptic and 32 controls placentas. All these samples had RIN values of >7. We used 40 ng RNA to make 48-plex RNA-seq libraries using modified STRT method with UMIs [48,49] and Globin lock method to deplete the abundant globin transcripts [50]. Briefly, RNA samples were placed in a 48-well plate in which a universal primer, template-switching oligos, and a well-specific 6 bp barcode sequence (for sample identification) were added to each well [51,52]. The synthesized cDNAs from the samples were then pooled into one library and amplified by single-primer PCR with the universal primer sequence. The libraries were sequenced with Illumina NextSeq 500 (RRID:SCR_014983) High Output (75 cycles).

The raw base call (BCL) files were demultiplexed using Picard tools (version 2.10.10; http://broadinstitute.github.io/picard/) ExtractIlluminaBarcodes and IlluminaBasecallsToSam to generate unaligned BAM files. These BAM files were converted to FASTQ files with Picard SamToFastq, and aligned to the human reference genome hg19, human ribosomal DNA unit (GenBank: U13369), and ERCC spike-ins (SRM 2374) using HISAT2 (version 2.1.0; RRID:SCR_015530) [53]. Aligned BAM files were then merged with the original unaligned BAM files to generate UMI-annotated BAM files using Picard MergeBamAlignment. These BAM files corresponding to each sample derived from four lanes were merged using Picard MergeSamFiles. Finally, potential PCR duplicates were marked by Picard Markduplicates. The resulting BAM files were processed with featureCounts (version 1.6.4) [54] to assign the reads to genes with options “-s 1 –largestOverlap –ignoreDup –primary”. Here, NCBI RefSeq (hg19) annotations downloaded from UCSC genome browser (RRID:SCR_005780) with addition of LINC02246 (chr21:16,189,159–16,290,989) were used. After removing three samples with extremely low number of mapped reads, the bias in read counts between the two libraries was corrected using the NBGLM-LBC method [55]. Differential expression analysis between 30 control and 32 severe preeclampsia samples was performed using DESeq2 (RRID:SCR_015687) as described above. The output ‘log2FoldChange’ of commonly detected 124 genes was used for the correlation with TAC-seq data. Because the FINNPEC patient consent doesn't allow sequencing data storage in public repositories, the STRT data are not publicly available.

qRT-PCR verification of placental IFNA1 expression. Because IFNA1 expression remained under the detection limit by STRT, we verified its expression by qRT-PCR. Each RNA sample (32 preeclamptic and 30 controls placentas; 800 ng per sample) used in STRT was reverse transcribed to cDNA in 20  μl reaction using the RT2 First Strand Kit (Qiagen) according to manufacturer's instructions. Real-Time PCR was performed using RT2 qPCR Primer Assays for IFNA1 (Qiagen; Cat# PPH01321B-200) and YWHAZ (Qiagen; Cat# PPH01017A-200), and RT2 SYBR Green Mastermixes. For the low expression of IFNA1, we used 5 µl of non-diluted cDNA template per reaction for IFNA1 and 2 µl for the housekeeping gene YWHAZ. Real-time PCR was performed in 25 µl reactions according to manufacturer's instructions using the Bio-Rad CFX96 qPCR Instrument (Bio-Rad Laboratories, Inc., Hercules, CA). The ΔΔCt method was used for data analysis and Welch's t-test to compare the relative expression of IFNA1 in preeclamptic (n = 29) versus control (n = 24) placentas with successful amplification.

We replicated the finding of IFNA1 upregulation in preeclamptic placentas using an independent publicly available dataset. RNA-seq data on expression profiles of 21 control and 20 preeclamptic placentas were downloaded from the GEO (RRID:SCR_005012), [42] under the accession number GSE114691 [56]. Normalization and differential expression analysis between these two groups were performed using DESeq2 (RRID:SCR_015687) with default parameters. Here, genes that had more than 10 read counts in total were assessed.

Immunohistochemistry and image analysis. We chose the HLA-G Monoclonal Antibody MEM-G/1 (Thermo Fisher Scientific; Cat# MA1–19,219; RRID: AB_1076732), recognizing both membrane-bound and soluble isoforms of HLA-G, for immunohistochemistry based on its lacking cross-reactivity with classical HLA I antigens [57]. The HLA-G antibody was diluted at 1:60 and rabbit monoclonal antibody to cytokeratin 7 (CK7) (Roche Diagnostics, Rotkreuz, Switzerland; Cat# 790–4462) at 1:10. To study cellular source of IFNA1 expression, mouse monoclonal antibody F-7 to interferon α (1/13) (Santa Cruz Biotechnology, Dallas, TX; Cat# sc-373,757; RRID:AB_10915749) was used and diluted at 1:100. Ultra-Sensitive ABC Peroxidase Staining Kit (Thermo Fisher Scientific Cat# 32052) was used according to the manufacturer's instructions followed by Metal Enhanced DAB Substrate Kit (Thermo Fisher Scientific; Cat# 34065).

The anti-HLA-G (n = 45) and anti-CK7 (n = 45) immunostained slides, in parallel with hematoxylin and eosin stained slides, were digitized with a whole-slide scanner (Pannoramic 250 FLASH, 3DHISTECH Ltd., Budapest, Hungary) equipped with a 20x objective (numerical aperture 0.8) and a 1x adapter, and a Global shutter CMOS camera with 4096 × 3072 pixels sized 5.5 μm x 5.5 μm (Adimec QUARTZ Q-12A180 camera) resulting in an image in which one pixel represents an area of 0.243 μm x 0.243 μm. Images were stored in a whole slide image format (MRX, 3DHISTECH Ltd., Budapest, Hungary) and further compressed to a wavelet file format (Enhanced Compressed Wavelet, ECW, ER Mapper, Intergraph, Atlanta, GA) with a target compression ratio of 1:10. The compressed virtual slides were uploaded to a whole-slide image management server (Aiforia Technologies Oy, Helsinki, Finland).

Two convolutional neural network-based algorithms were trained using a cloud-based software (Aiforia Create, Aiforia Technologies Oy, Helsinki, Finland) for image analysis to separately identify anti-CK7 and anti-HLA-G immunostained tissue areas. Each algorithm consists of two algorithms in sequence, the first segmenting the extravillous tissue areas and the second algorithm to identify trophoblasts that exhibit a strong HLA-G expression or a strong CK7 expression respectively. The network for the detection of anti-HLA-G stained tissue areas was trained with 588 image annotations corresponding to an area of 0.0066 mm2 and 0.001 mm2 extravillous tissue. The algorithm for detecting anti-CK7 immunostained tissue areas was trained with 193 image annotations corresponding to an area of 0.11 mm2 and 22 mm2 extravillous tissue. We used a feature size of 15 μm for training the areas with strong HLA-G expression and a feature size of 10 μm for training the algorithm for detection of anti-CK7 stained tissue areas. The augmentation parameters for detection of anti-HLA-G immunostaining were 0°–360° rotation, ±20% aspect ratio change, ±10% shear distortion, ±20% brightness and ±20% contrast change and augmentation parameters for the detection of CK7 expression were 0°–360° rotation, ±5% aspect ratio change, ±50% shear distortion, ±5% brightness and ±5% contrast change. All data were flipped both vertically and horizontally. The algorithm performances were validated against manual annotations done by one of the researchers (NL).

HLA-G/CK7 protein ratios were received from the image analysis. For the placentas with two biopsies, mean ratios were calculated and used in the analyses. HLA-G/CK7 protein ratios were plotted against HLA-G/CK7 mRNA ratios, which were calculated for each placenta from the TAC-seq data, after adjustment the expression data by offspring sex, mode of delivery, and gestational age at birth (<32, 32–36, ≥37 weeks) as described. Spearman correlation coefficient was calculated to study the relationship between HLA-G/CK7 mRNA and protein expression. HLA-G/CK7 protein ratios were studied in relation to the gestational weeks at delivery. Comparison of nonlinear fits in preeclamptic and controls pregnancies was performed. We further studied the differences in protein expression between the two sexes in preeclampsia and control offspring by using the Student's t-test.

ELISA. Human IFNα ELISA Kit (Sigma-Aldrich, Merck, Darmstadt, Germany; Cat# RAB0541) was used according to the manufacturer's instruction to detect IFNα−1/13 in third trimester maternal serum samples (n = 209; 158 preeclampsia and 51 controls), taken at variable time points during hospitalization (0 to 31 days before the delivery), and fetal cord plasma samples (n = 176; 66 preeclampsia and 110 controls), taken immediately after delivery. We first tested the mean absorbance of two duplicates from a total of 32 samples, and found the average coefficient of variance (cv) of 6%, which was below the intra-assay cv of <12% reported by the manufacturer. For the limited availability of samples, the final experiments were run without duplicates, using 100 µl of each sample, SYNERGY H1 Microplate reader, and the Gen5 Software (BioTek Instruments, Winooski, VT). Standard curve was defined according to the kit's instructions and standard dilutions. Values were converted to ng/mL by reference to a standard curve which was generated in parallel with the test samples in each run.

Maternal and fetal IFNα concentrations were compared by the Mann-Whitney U test, and studied in relation to the disease status by the Chi-Square test in preeclamptic vs. non-preeclamptic pregnancies. Levels of maternal blood pressure were compared between the subgroups with IFNα concentrations below and at/above the median value (Mann-Whitney U test).

Statistics. Statistical analyses were carried out using GraphPad Prism software, version 8.3.0 (GraphPad Software, San Diego, CA), unless otherwise stated. Two-sided P values of 0.05 or less were considered statistically significant.

Chi-Square test was used to compare the proportions of males and females at birth in preeclamptic and control pregnancies, with and without history of miscarriages, and Spearman correlation coefficients to study the relationship between the sex ratio and gestational ages at birth.

The sex ratios were further analysed using a generalized linear model (GLM; events-trials, binomial distribution and logit link function) in SPSS statistical software (version 22; IBM, Armonk, NY). For the 1000 Genomes dataset, the model included the allele 1, allele 2, location (African vs. Non-African populations) and haplotype frequency as main effects, and the interaction of allele 1 and allele 2, as well as the interaction of location and haplotype frequency. For the FINNPEC dataset, the final model included the allele 1, allele 2, haplotype frequency and maternal condition (preeclampsia or control) as main effects, and the interaction of allele 1 and allele 2. In these models, the dependent variable (sex ratio) was analysed as the number of males out of total offspring.

The numbers of homozygous and heterozygous individuals for HLA-G 3′UTR and the associated ORs and 95% CIs were calculated by Chi-Square test. The numbers of the observed HLA-G 3′UTR haplotypes were compared between stillborn fetuses and controls, and in preeclampsia offspring and controls, using Chi-Square test and comparing the haplotypes one by one against all other haplotypes.

RNA and protein expression analyses and statistical methods are detailed in the corresponding sections describing the methods for TAC-seq, STRT RNA-seq, RT-qPCR, immunohistochemistry, and ELISA.

Ethics Statement and Study approval. Written informed consent was received from all participants in the FINNPEC cohort [26]. Samples for the stillbirth cohort were from the Biobank Karolinska. All methods were carried out in accordance with protocols approved by the participating institutions (Helsinki University Central Hospital 149/E0/07; Karolinska Institutet 2008/670–31/2 and 2015/1704–31/1) and in agreement with the Declaration of Helsinki.

Results

Male fetuses are underrepresented in preeclamptic births

We examined birth sex ratios (male/female) in a Finnish population cohort of 1.79 million live and stillbirths. This cohort included 38,752 preeclamptic births (2.2%). The sex ratio was 1.02 in preeclamptic and 1.05 in non-preeclamptic births (OR 0.97, 95% CI: 0.95 to 0.99; P = 0.006; Chi-Square test), and the earlier the gestational age at birth, the lower the number of male offspring in preeclampsia (Spearman correlation coefficient, 0.80; P = 0.007) (Fig. 1a). The sex bias was in striking agreement with Norwegian data on 1.82 million births and 44,000 preeclamptic pregnancies, [16] and with others [18]. Together with the data on unbiased human sex ratio at conception [1] and familial association of preeclampsia and miscarriages, [58] our results suggested male fetuses being lost before the onset of early preeclampsia, or linked to a higher risk of term disease, or both.

Fig. 1

Biased birth-sex ratio in preeclampsia is associated with balancing selection of3′UTR diplotypes. (a) Population-level birth sex ratios in preeclampsia and control pregnancies. Preeclamptic births after 41 weeks were rare (n = 512) and the sex ratio of 0.99 comparable with that of the controls (data not shown). (b) The number of previous pregnancies and associated birth sex ratios in the population, in the cohort of preeclamptic pregnancies, and in the subgroup of preeclamptic pregnancies with small for gestational age (SGA) offspring. Numbers of individuals are represented in the bars and P values (Chi-Square test) for birth sex differences above the graphs (a and b); dashed lines indicate the population birth sex ratio of 1.05 and n.s. denotes non-significant. The distribution of HLA-G 3′UTR diplotypes (e.g. 1/2 denoting UTR-1/UTR-2) and the associated birth sex ratios in (c) preeclampsia, (d) stillborn fetuses, and (e) control offspring. Horizontal lines indicate the balanced sex ratio of 1.0. The preeclampsia cohort (c and e) indicates that the HLA-G diplotype is the main determinant of the bias in the birth sex ratio (GLM; allele 1: F5,21=0.44, P = 0.82; allele 2: F7,21=3.51, P = 0.01; allele 1 by allele 2 interaction: F13,21=9.15, P<0.001; haplotype frequency: F1,21=0.67, P = 0.42; maternal condition: F1,21=3.58, P = 0.07). (f) The 1000 Genomes dataset indicates that the haplotype, haplotype frequency and geographic location all explain the bias in the sex ratio between (g) African and (h) Non-African individuals (GLM; allele 1: F7,24=41.63, P<0.001; allele 2: F8,24=8.23, P<0.001; allele 1 by allele 2 interaction: F21,24=2.27, P = 0.028; location: F1,24=6.28, P = 0.02; haplotype frequency: F1,24=7.87, P = 0.01; location by haplotype frequency interaction: F1,24=8.47, P = 0.008).

In the preeclampsia cohort (n = 1,538), the higher number of previous pregnancies was associated with reduced number of male offspring in preeclamptic women, especially in those delivered small for gestational age (Fig. 1b). This deficiency of male offspring was not present in first pregnancies complicated with preeclampsia. The deficiency of males in multiparous preeclamptic births could have two explanations: (1) the 'missing males' are in the non-preeclamptic birth category (maleness is protective) or (2) males in potentially preeclamptic pregnancies are lost before birth. To provide evidence of preferential loss of males, we studied population-level data on women with one or more miscarriages and found the reduced sex ratio of 0.98 in consecutive preeclamptic (n = 6,469) compared with 1.04 in non-preeclamptic births (n = 355,683) (OR 0.95, 95% CI: 0.90–0.99; P = 0.03; Chi-Square test). Moreover, we studied a Swedish stillbirth cohort (n = 277), which showed increased intrauterine loss of male fetuses with the sex ratio of 1.18, being in agreement with earlier findings. [1] Altogether, these findings of sex-ratio distortion propose that maternal immune responses to the male-specific Y chromosomal histocompatibility (H-Y) antigen, [59] present on fetal cells or microchimeric cells retained in the mother's body from the previous pregnancy, [60] contribute to fetal survival and preeclampsia.

Fetal HLA-G 3′UTR modulates human sex ratio

HLA-G downregulation is linked to its 3′ untranslated regulatory region (3′UTR) (Fig. S1), which is associated with distinct full-length haplotypes, and modulates mRNA stability and decay, microRNA targeting, and splicing [5,22,61,62]. To uncover the link between HLA-G and the sex ratio, we studied 3′UTR haplotype pairs, diplotypes, representing functional units. Rare diplotypes showed opposite patterns of sex ratios in preeclampsia offspring (Fig. 1c) and stillborn fetuses (Fig. 1d), and HLA-G diplotype was the main determinant of the birth sex bias between preeclampsia and control offspring (Fig. 1e). Similarly, HLA-G haplotypes were associated with the bias in the sex ratio between African versus Non-African populations in the 1000 Genomes data (Fig. 1f, 1g and 1h), the African populations having the highest global prevalence of preeclampsia [63]. Altogether, rare diplotypes showed more variation in the sex ratio than the common ones, which were close to the balance of 1.0. These results indicate balancing selection and negative frequency-dependence of allelic fitness (rare allele advantage), meaning that rare alleles are favored and the fitness of each allele decreases as its frequency increases [64].

We further tested segregation of individuals homozygous for HLA-G 3′UTR. Decreased HLA-G homozygosity was observed in African (18.3% homozygous) versus Non-African populations (24.4%), and in the stillbirth cohort (18.4%) versus population controls (25.7%), without sex-specific differences. Preeclampsia (26.8%) versus control offspring (26.3%) showed no differences in the overall homozygosity. However, homozygous males (31.7%; 35.2% in first pregnancies) versus females (22.3%; 21.4% in first pregnancies) were overrepresented in preeclamptic births (Fig. 2). These findings provide more support for balancing selection and are suggestive of heterozygote advantage. This is a special case of negative frequency-dependent selection and means that a greater proportion of the rarer allele than the more common allele is present in advantaged heterozygotes [2,65].

Fig. 2

Patterns of3′UTR homozygosity support balancing selection. The odds ratios (dots) and 95% confidence intervals (whiskers) for homozygous versus heterozygous individuals for HLA-G 3′UTR in African versus Non-African populations, stillborn fetuses versus populations controls, and in preeclampsia versus control offspring. The comparisons between homozygous and heterozygous males versus females were performed in African individuals, stillborn fetuses, and in preeclampsia offspring. Only P values <0.05 are shown (Chi-Square test).

HLA-G 3′UTR shows association with preeclampsia and stillbirth

The two major HLA-G 3′UTR haplotypes, UTR-1 and UTR-2, differ for five polymorphic sites, mostly associate with the same HLA-G protein, and have global frequencies of 20–30% each [22]. Of them, the evolutionarily most recent UTR-1 showed reduced frequency (P = 0.01; OR 0.77, 95% CI: 0.63–0.94; Chi-Square test) and the ancestral haplotype UTR-2 [66] increased frequency (P = 0.01; OR 1.34, 95% CI: 1.07–1.69; Chi-Square test) in preeclampsia versus control offspring in first-time mothers. For UTR-2, the association remained when offspring from later pregnancies were included (P = 0.04; OR 1.22, 95% CI: 1.01–1.47; Chi-Square test) (Table S3). The observed ORs are considered good approximations of the risk ratio, because of the low preeclampsia rate of 2.1% during the preeclampsia series collection. The opposite effects of UTR-1 and UTR-2 on preeclampsia risk in first pregnancies support the advantage of divergent alleles and heterozygosity. Moreover, previous studies link UTR-2 to low expression, immune-mediated disorders, and pregnancy complications, and UTR-1 to high HLA-G expression and increased mRNA stability [5,11,67]. In line with this, the stillbirth cohort showed a similar but more significant protective effect of UTR-1 (P = 0.0001; OR 0.65, 95% CI: 0.53–0.81; Chi-Square test). Conversely, UTR-5 was detected as a risk haplotype for stillbirth (P = 0.03; OR 1.72, 95% CI: 1.06–2.79; Chi-Square test) (Table S4). UTR-5 is considered the oldest of all haplotypes and is closest to orthologous sequences in non-human primates [66]. We further confirmed that the 14 bp insertion polymorphism, previously implicated in pregnancy complications and present in haplotypes UTR-2, −5 and −7, [11,22] showed a modest association with both stillbirths (P = 0.04; OR 1.24, 95% CI: 1.01–1.52) and preeclampsia (P = 0.03; OR 1.21, 95% CI: 1.02 to 1.44; Chi-Square test).

HLA-G downregulation and IFNA1 upregulation characterize preeclamptic placentas

To test the hypothesis of maternal immunoreactivity as the mechanism of fetal selection, we studied placental samples (n = 163) for mRNA expression of 136 HLA-G related, placenta-specific, and preeclampsia-associated genes (Table S1) using the highly quantitative TAC-seq method [41] and validated by STRT RNA sequencing (Fig. S3). HLA-G expression, measured by 15 different probe pairs, showed consistent downregulation in preeclamptic placentas. HLA-G 3′UTR diplotypes were associated with differential HLA-G expression in controls but HLA-G expression was always low in preeclamptic placentas, irrespective of the haplotype combinations (Fig. 3a and Fig. S2). Although UTR-1/UTR-1 was among those with the highest and UTR-2/UTR-2 among those with the lowest HLA-G expression in non-preeclamptic placentas, differences between the individual haplotypes were non-significant due to the small number of samples. To support the association between HLA-G 3′UTR, preeclampsia, and fetal sex, the diplotypes producing more males were associated with higher placental HLA-G mRNA expression (Fig. 3b and 3c). This finding suggested that sex-related segregation of HLA-G haplotypes is associated with differential HLA-G expression levels. Similarly, we observed higher HLA-G protein expression in preeclamptic placentas from male compared with female offspring (Fig. 3d). We further confirmed the correlation between HLA-G/CK7 mRNA and protein expression (Fig. 3e and 3f). Collectively, these findings supported sex-related differences and downregulation of HLA-G in preeclamptic placentas independently of trophoblast count, as defined by CK7 expression.

Fig. 3

downregulation in preeclamptic placentas. (a) HLA-G mRNA levels and 3′UTR diplotypes in preeclamptic and control placentas (n = 163; P = 1.63 × 10−4; Welch's t-test); comparisons between diplotypes remained non-significant. The number of samples are shown above the graph and linear regression lines separately for controls (blue) and preeclampsia (red). Haplotypes are arranged based on the highest to the lowest expression observed in the control placentas. (b) Linear regression summarizes HLA-G diplotypes with positive residuals, defined as male excess, and with negative residuals, defined as female excess, in the cohort of preeclamptic and control offspring (n = 1,247). Analysed based on this classification, (c) demonstrates higher placental (n = 163) HLA-G mRNA expression in association with diplotypes showing male excess (Welch's t-test). (d) HLA-G/CK7 protein ratios, and two regression lines with significantly different slopes (P = 0.02), demonstrate HLA-G protein expression in preeclamptic (n = 13; red) and control placentas (n = 11; blue). HLA-G protein expression was lower in female (n = 7) compared with male offspring (n = 6) in preeclampsia (PE) (Student's t-test), but no sex-related differences arose in controls. Horizontal line indicates the mean of all samples. (e) Spearman correlation (r) visualized as scatterplot shows HLA-G/CK7 protein and mRNA expression in 24 placentas. RNA expression levels are shown adjusted by offspring sex, mode of delivery and gestational age at birth; only panel c shows unadjusted data. (f) HLA-G immunoreactivity (marked in red by the image analysis) on extravillous trophoblasts in the representative preeclamptic (PE) placenta (female born at 29 weeks; HLA-G UTR-2/UTR-4) with 2% and the control (Ctrl) placenta (female born at 32 weeks; HLA-G UTR-1/UTR-2) with 13% HLA-G positive trophoblasts (scale bars, 500 µm). For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.))

Altogether 36 genes were differentially expressed between severe preeclampsia (n = 81) and controls (n = 63); 8 upregulated and 28 downregulated (Fig. 4a; Table S5). Collectively, we found downregulation of HLA-G and other markers for differentiated trophoblasts (IL2RB, MMP2, SPRR2G, LAIR2, TGFB2), [29] HLA-G receptors (LILRB1 and LILRB2), [5] immune escape molecules of the B7 family (PDCD1LG2/B7-DC/CD273 and CD276/B7-H3), [32] Th2 cytokine IL10, [68] and many leukocyte markers (CD300E, MS4A7, AQP9, GNLY) in preeclamptic placentas. We found no evidence that other immune-related genes such as IDO, HLA-E or HLA-F would compensate HLA-G downregulation in preeclampsia. Altogether, we observed an overall downregulation of tolerogenic genes in preeclamptic compared with control placentas.

Fig. 4

Downregulation of tolerogenic genes and upregulation ofin preeclamptic placentas. (a) Volcano plot of downregulated (blue dots) and upregulated (red dots) placental genes (Padj<0.05; DESeq2 Wald test) in severe preeclampsia (n = 81) versus controls (n = 63). Horizontal dashed line indicates Padj=0.05 and vertical dotted line indicates no fold change. (b) ROC curves show a significantly higher mean (±SE) area under the curve for HLA-G and IFNA1 as predictors of preeclampsia when compared with FLT1 and PGF. Expression levels were adjusted by offspring sex, mode of delivery, and gestational age at birth. For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.

Fms-like tyrosine kinase 1 (FLT1), a major driver of maternal hypertension in preeclampsia, [69] showed the second highest upregulation in preeclamptic placentas. Notably, a major IFN subtype alpha-1 (IFNA1), [70] a cytokine modulating innate and adaptive immune responses, [71] was the most upregulated gene in preeclampsia (Fig. 4a), while interferon beta and gamma showed no differences. We further confirmed upregulation of IFNA1 in preeclampsia by qRT-PCR (Fig. S4) and in an independent publicly available RNA sequencing dataset (Fig. S5) [56]. We studied cellular source of IFNA1 and detected protein expression in fetal trophoblasts (Fig. S6). We further observed that fetal, but not maternal, circulating IFNα levels were associated with maternal hypertension (Fig. S7). These findings are in line with the absence of transplacental transfer of IFNα [72] and support the fetal compartment in IFNα production.

In our TAC-seq dataset, placental expression of HLA-G and IFNA1 distinguished preeclamptic placentas better than expression of the current biomarkers FLT1 and placental growth factor (PGF) (Fig. 4b) [73]. We further studied gene expression patterns related to the HLA-G regulatory haplotypes, and found that, FLT1 was the main measure that differed between common versus rare, and homozygous versus heterozygous, placentas with respect to control and preeclampsia (Fig. 5).

Fig. 5

3′UTR haplotypes show asymmetry forexpression.Shown are expression levels of HLA-G, IFNA1 and FLT1 divided based on the HLA-G 3′UTR haplotypes targeted by selection. (a) The samples (n = 95) with the most common diplotypes (UTR-1/UTR-4, UTR-1/UTR-2, UTR-1/UTR-1, and UTR-2/UTR-4) and the samples (n = 68) with rare diplotypes (all others). (b) The heterozygous (n = 120) samples and the homozygous samples (n = 43) for HLA-G 3′UTR. Expression levels were adjusted by offspring sex, mode of delivery, and gestational age at birth. P values were calculated by the Welch's t-test.

Discussion

Our study provides evidence that failure in fetal HLA-G mediated immune escape is associated with stillbirth and preeclampsia, and that intriguingly, the birth sex ratio in these pregnancy complications reflects the action of this mechanism in vivo. We link regulatory HLA-G haplotypes to the offspring sex and pregnancy complications by the two major mechanisms of balancing selection at the MHC locus: negative frequency-dependent selection (rare allele advantage) and heterozygote advantage [64,65]. While negative frequency-dependent selection has been proposed to keep sex ratios balanced, [74] there is very limited earlier evidence for a genetic link to sex ratio. Thus, to our knowledge, we present the first evidence in any species of a gene modulating the sex ratio and suggest that fetal HLA-G is a sexually antagonistic locus determining maternal-fetal interaction and human pregnancy success. Furthermore, we provide evidence that balancing selection is acting on contemporary human populations and contributes to human pregnancy complications. Moreover, our findings add significant novel data to support the hypothesized role of aberrant IFN signaling in pregnancy complications, [75] as shown in mouse Zika virus infection, [76] but with limited previous data on human pregnancies.

Although we studied large cohorts, the heterogeneity of preeclampsia and challenges in the establishment of explicit diagnostic criteria remain a matter of speculation. Similarly, the 2.2% rate of preeclampsia in the Finnish birth registry data is close to the lowest global estimate of 2%, [13] but due to the general and equal access to maternal care with globally low infant mortality may miss only a few cases. Our registry data involved only births after 22 weeks of gestation, thus providing only indirect evidence of the association between male fetal loss and the biased birth sex ratio in preeclampsia. However, our results are broadly compatible with the earlier findings of female-biased birth sex ratio in early preeclampsia [16], [17], [18] and with greater losses of male conceptions before female live births with preeclampsia or low birth weight, as demonstrated in couples with recurrent miscarriages [77,78]. Furthermore, both clinical cohorts supported the association between HLA-G haplotypes, the sex ratio, and pregnancy complications. That even the unselected global 1000 Genomes populations showed association between HLA-G haplotype segregation and the sex ratio, further supports the robustness of these results.

Our results show underrepresentation of male offspring in preeclamptic births, especially those delivered preterm or small for gestational age. This finding supports the vulnerability of male fetuses to maternal inflammation [79] and late miscarriages and stillbirths [1]. Because there was no sex-ratio distortion in first pregnancies, but the deficiency of males was present only in multiparous births, our findings may reflect a different etiology of preeclampsia in first and subsequent pregnancies. Our results show, however, the highest risk of preeclampsia in first pregnancies with a male fetus homozygous for HLA-G 3′UTR, thus supporting the association between male fetal sex and preeclampsia also in first-time mothers. Furthermore, our findings indicate that fetal HLA-G might be the locus mediating a deficiency of HLA homozygotes through maternal-fetal interaction [2] and explaining HLA similarity in couples with recurrent miscarriages [3]. In contrast, heterozygote advantage looks apparent for African populations, but heterozygosity also shows an association with stillbirth in the Swedish series. This suggests that heterozygote advantage might be related to special circumstances, such as resistance to infectious diseases as proposed for the MHC locus [65]. Thereby, balancing selection of fetal HLA-G alleles points toward pathogen interactions and the suggested association between viral infections and preeclampsia [80,81].

Our results strengthen the earlier findings of HLA-G downregulation in preeclamptic placentas [6], [7], [8]. We studied a large number of samples and HLA-G/CK7 ratios to demonstrate that HLA-G mRNA and protein downregulation in the preeclamptic placentas was not explained by the lower number of trophoblasts in the samples. Although our dataset was small to show significant differences between specific HLA-G 3′UTR haplotypes and expression levels, we found an overall trend and added novel data by demonstrating the association between HLA-G mRNA and protein expression and fetal sex. Additionally, we showed downregulation of HLA-G receptors, such as LILRB1 linked to lower expression in first pregnancies, [82] and many other tolerogenic genes in preeclamptic placentas. Furthermore, our results suggest that IFNα signaling might drive loss of fetal immunotolerance in preeclampsia, as implicated in systemic lupus erythematosus (SLE), [83] in tumor rejection, [84] and in solid organ transplantation [85]. To support this, increased IFNα/β signaling is associated with preeclampsia in SLE women, [86] inflammation-driven preterm birth in mice and humans, [87,88] and in a pathway crosstalk analysis of preeclampsia [89]. While viruses and bacteria are the main triggers of IFNα/β production but rarely present in the human placenta, [90,91] even fetal DNA immune complexes alone, excessively released from damaged cells of the preeclamptic placenta, [92,93] may activate the autophagy pathway and promote prolonged IFNα production [83,94]. Importantly, hydroxychloroquine (HCQ) reduces IFNα production and preeclampsia in SLE women, [95], [96], [97] improves endothelial function in an in vitro model of preeclampsia [98] and in SLE mice, [99] as well as prevents placental inflammation and fetal loss in mice [100,101]. Therefore, the mechanisms presented here encourage targeting IFNα signaling in clinical trials and testing whether HCQ has general therapeutic effects in preeclampsia.

This study indicates that HLA-G heterozygosity, as a determinant of fetal survival, might be associated with enhanced FLT1 responses. Remarkably, previous studies link overexpression of FLT1 to fetal self-defense against decidual VEGF signaling [102]. Supporting fetal defense and maternal-fetal conflict, high fetal FLT1 levels confer resistance against placental malaria and fetal loss [103] and maternal (decidual) macrophages express VEGF in placental malaria [104]. Considerably, maternal hypertension and elevated sFLT1 levels are not specific to preeclampsia but arise in first-time mothers with placental malaria [104]. Moreover, low HLA-G expression emerges in preeclamptic placentas, but also in trophoblasts infected by Plasmodium falciparum [105] and some viruses (CMV, HSV) [106], [107], [108]. In further support of shared mechanisms, placental malaria and preeclampsia show high rates and seasonal co-occurrence in Africa,[109] where both the HLA-G low-expressing “null allele” and high-expressing FLT1 alleles are under positive selection, [103,110] and the frequency and sex ratios associated with HLA-G 3′UTR haplotypes differ from other populations as shown in this study. Furthermore, the African “null allele” is mostly associated with UTR-2, [22] the haplotype showing an increased risk of preeclampsia in this study. Thus, we suggest that low HLA-G expression may be advantageous in pathogen-rich environments to mount protective maternal immune responses, while concurrently detrimental, causing unfavorable maternal anti-fetal immune reactions, fetal loss, and preeclampsia. In agreement with this scenario, ancestral fetal HLA-G haplotypes showed associations with preeclampsia and stillbirths in our series, thereby supporting their disadvantage in modern environments. Furthermore, high HLA-G expression is associated with higher implantation rates [111] and uneventful pregnancies as presented here, but might confer an increased risk of offspring to malaria and other infections [112,113].

Collectively, our results suggest that balancing selection at fetal HLA-G locus modulates fetal survival, preeclampsia and birth sex ratio. We propose that preeclampsia involves coevolutionarily balanced mechanisms: downregulation of fetal HLA-G mediated immunotolerance, upregulation of maternal immune responses as reflected by male fetal loss, and upregulation of fetal self defense mechanisms such as overexpression of FLT1 and IFNA1. Although these mechanisms may have evolved to protect both the fetus and the mother from invading pathogens, they increase the risk of maternal immune attack against the fetus, fetal loss, and preeclampsia. We suggest that this coevolutionary balancing explains the highest incidence of preeclampsia in sub-Saharan Africa and promotes genetic diversity in fetal HLA-G. Future studies should test whether other mechanisms, such as the lack of alleles in the fetus but their presence in the mother, [114,115] might modulate fetal selection and contribute to human pregnancy complications. While HLA-G shows physiological expression almost uniquely in trophoblasts, its high expression that benefits the viability of fetuses during pregnancy may trade-off with HLA-G neoexpression, promoting immune evasion of viruses and malaria as well as tumors later in life [5,116]. Therefore, the benefit of differential HLA-G expression may be dependent on the pathogen load and the tissue of expression, providing a role for balancing selection on immune-related traits beyond pregnancy in human evolution.

Declaration of Competing Interest

Dr. Katayama reports grants from Jane and Aatos Erkko Foundation, grants from Swedish Research Council, during the conduct of the study. Dr. Laivuori reports grants from Jane and Aatos Erkko Foundation, grants from Päivikki and Sakari Sohlberg Foundation, during the conduct of the study; grants from The Competitive State Research Financing of the Expert Responsibility area of Tampere University Hospital, grants from Tampere University Hospital Foundation, outside the submitted work. Drs. Krjutškov and Kere have a patent on the TAC-seq method pending. Other authors have nothing to disclose.

Author contributions

SW and JK conceived and designed the study. ES, EI, KP, EK, SH, HL and JK collected the patient cohorts and prepared the samples. SW, HS, IMH, SE, TS and KK performed the experiments. SW, MY, HT, MG, SK, JW, NL, JL, MM and JK generated, analysed, and interpreted the data. SW, MY and JK wrote the manuscript and all authors substantially revised the manuscript.