# Pan-Genomic Regulation of Gene Expression in Normal and Pathological Human Placentas

Clara Apicella, Camino Ruano, Basky Thilaganathan, Asma Khalil, Veronica Giorgione, Géraldine Gascoin, Louis Marcellin, Cassandra Gaspar, Sébastien Jacques, Colin Murdoch, et al.

► To cite this version:

Clara Apicella, Camino Ruano, Basky Thilaganathan, Asma Khalil, Veronica Giorgione, et al. Pan-Genomic Regulation of Gene Expression in Normal and Pathological Human Placentas. Cells, 2023, 12 (4), pp.578. 10.3390/cells12040578. hal-04096063

HAL Id: hal-04096063

https://u-paris.hal.science/hal-04096063

Submitted on 23 Oct 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Article

# Pan-Genomic Regulation of Gene Expression in Normal and Pathological Human Placentas

Clara Apicella<fnref n="1" />, Camino S. M. Ruano<fnref n="1" />, Basky Thilaganathan <fnref n="2" />,<fnref n="3" /> ID, Asma Khalil <fnref n="2" />,<fnref n="3" /> ID, Veronica Giorgione <fnref n="2" />,<fnref n="3" /> ID,
Géraldine Gascoin <fnref n="4" />, Louis Marcellin <fnref n="5" /> ID, Cassandra Gaspar <fnref n="6" />, Sébastien Jacques <fnref n="1" />, Colin E. Murdoch <fnref n="7" /> ID,
Francisco Miralles <fnref n="1" />, Céline Méhats <fnref n="1" /> ID and Daniel Vaiman <fnref n="1" />,* ID

<figure><img src="image_2.png" /></figure>

check for
updates

**Citation:** Apicella, C.; Ruano, C.S.M.; Thilaganathan, B.; Khalil, A.; Giorgione, V.; Gascoin, G.; Marcellin, L.; Gaspar, C.; Jacques, S.; Murdoch, C.E.; et al. Pan-Genomic Regulation of Gene Expression in Normal and Pathological Human Placentas. *Cells* **2023**, *12*, 578. https://doi.org/10.3390/cells12040578

Academic Editors:
Hironori Takahashi, Osamu Samura
and Koumei Shirasuna

Received: 4 November 2022

Revised: 17 January 2023

Accepted: 28 January 2023

Published: 10 February 2023

<figure><img src="image_3.png" /></figure>

**Copyright:** © 2023 by the authors.
Licensee MDPI, Basel, Switzerland.
This article is an open access article
distributed under the terms and
conditions of the Creative Commons
Attribution (CC BY) license (https://
creativecommons.org/licenses/by/
4.0/).

**Abstract:** In this study, we attempted to find genetic variants affecting gene expression (eQTL = expression Quantitative Trait Loci) in the human placenta in normal and pathological situations. The analysis of gene expression in placental diseases (Pre-eclampsia and Intra-Uterine Growth Restriction) is hindered by the fact that diseased placental tissue samples are generally taken at earlier gestations compared to control samples. The difference in gestational age is considered a major confounding factor in the transcriptome regulation of the placenta. To alleviate this significant problem, we propose here a novel approach to pinpoint disease-specific cis-eQTLs. By statistical correction for gestational age at sampling as well as other confounding/surrogate variables systematically searched and identified, we found 43 e-genes for which proximal SNPs influence expression level. Then, we performed the analysis again, removing the disease status from the covariates, and we identified 54 e-genes, 16 of which are identified de novo and, thus, possibly related to placental disease. We found a highly significant overlap with previous studies for the list of 43 e-genes, validating our methodology and findings. Among the 16 disease-specific e-genes, several are intrinsic to trophoblast biology and, therefore, constitute novel targets of interest to better characterize placental pathology and its varied clinical consequences. The approach that we used may also be applied to the study of other human diseases where confounding factors have hampered a better understanding of the pathology.

**Keywords:** placenta; preeclampsia; expression Quantitative Trait Loci

## 1. Introduction

Eutherian mammals made a drastic evolutionary 'choice' for survival, the one to host the embryo and to develop the fetus over a relatively extended period of time. Since the feto-placental unit is, in fact, a semi-allograft, that in all other circumstances should be rejected by the maternal immune system. This tolerance period ranges from a few days in marsupials (13 days in the Virginia opossum, where poor coping with inflammation apparently shortens the length of the feto-maternal tolerance period [1]) to 22 months in the African elephant. The pivotal organ for this heavy evolutionary option is the placenta,

<footnote n="1">¹ Team 'From Gametes to Birth', Institut Cochin, U1016 INSERM, UMR 8104 CNRS, Paris-Descartes University, 75014 Paris, France</footnote>

<footnote n="2">² Fetal Medicine Unit, St George's University Hospitals NHS Foundation Trust, London SW17 0RE, UK</footnote>

<footnote n="3">³ Vascular Biology Research Centre, Molecular and Clinical Sciences Research Institute, St George's University of London, London SW17 0RE, UK</footnote>

<footnote n="4">⁴ Department of Neonatology, Angers University Hospital, F-49000 Angers, France</footnote>

<footnote n="5">5 Department of Gynaecology, Obstetrics and Reproductive Medicine, Centre Hospitalier Universitaire (CHU) Cochin Faculté de Médecine, Assistance Publique-Hôpitaux de Paris (AP-HP), Hôpitaux Universitaires Paris Centre (HUPC), Université de Paris, 138 Boulevard de Port-Royal, 75014 Paris, France</footnote>

<footnote n="6">⁶ Sorbonne Université, Inserm, UMS Production et Analyse des données en Sciences de la vie et en Santé, PASS, Plateforme Post-génomique de la Pitié-Salpêtrière, 75013 Paris, France</footnote>

<footnote n="7">⁷ Systems Medicine, School of Medicine, University of Dundee, Dundee DD1 9SY, UK</footnote>

<footnote>* Correspondence: daniel.vaiman@inserm.fr; Tel.: +33-1-44412301; Fax: +33-1-44412302</footnote>
which manages immunotolerance, exchange of nutrients, withdrawal of waste, and pro-
gestational hormone production. The inherent complexity associated with this role of the
interface makes its regular function a fragile equilibrium, which in humans evolve relatively
often towards placental diseases, such as Hypertensive Disorders of Pregnancy (HDP),
Intra-uterine Growth restriction (IUGR), Gestational diabetes or other diseases related to
abnormal implantation [2,3]. Despite this obvious fundamental role, understanding the
genetic function of the human placenta is a difficult challenge, and overall this organ is not
explored as much as other organs. For instance, the Genotype Tissue Expression (GTex)
project, for instance, provides a comprehensive list of expression-Quantitative Trait Loci
(eQTLs), as well as splicing-QTLs (sQTLs) in 49 human tissues, that do not include the
placenta [4] Recent progress using global approaches made it possible to better analyze the
regulation of gene expression in the pathological placenta, such as in preeclampsia and
Intra-Uterine Growth Restriction [5,6]. Genetic regulation of gene expression by genetic
variants has been systematically analyzed in three recent articles that initiated deciphering
the landscape of genetically controlled genes in the placenta [7–9]. Our data presented
herein revisit the question of eQTL in the placentas and uses a novel methodological
approach to identify actual differences between normal and pathological placentas often
blurred by a major confounding factor, such as placental age. For this, we gave a strong
weight to the use and the definition of covariates. Covariates were defined by available
clinical (i.e., age of the placenta, disease status) and technical parameters (i.e., sample
preparation, origins), by principal components discovered from the transcriptome analysis,
and through the systematic identification of Surrogate Variables [10].

Since confounding factors are a major issue in placental genetics (especially the gesta-
tional age, the mode of delivery–C-section or Natural), we develop a novel approach that
could be used to identify relevant disease-associated genetic variants (SNPs) associated
with alterations of gene expression. In the first step, the QTL analysis was performed con-
sidering the disease status as a covariate, and then, in the second step, we did not correct
for disease status looking for eQTLs that might be influenced by genetics and disease and
being, therefore, potentially interesting in the context of placental pathologies. We surmise
that the additional couple of SNP-gene found were mostly and uniquely specific to the
disease. While it is not possible to affirm that the gene expression alterations are causal
to the disease, they could later be evaluated in other human cohorts of placental DNAs to
assess their potential predictive value heuristically.

The approach is summarized in Figure 1. This alternative analysis strategy could be
applicable to a wide range of eQTL studies, especially those in which gene expression is
believed to be highly influenced by clinical and technical variables and that is limited by
small sample sizes, such as those including disease samples.
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Overview of the experimental design. SNP = Single Nucleotide Polymorphism; PCA = Principal Component Analysis; SVA = Surrogate Variable Analysis; eQTL = expression Quantitative Trait Locus.</figcaption></figure>

**2. Methods**

**2.1. Summary of the Principle**

In the first step, we performed the analysis including the disease status (Control, Pre-eclampsia (PE), IUGR, and PE + IUGR) as a covariate, which allowed us to identify eQTLs that influence placental gene expression independently of the disease. In order to restrict the number of input genes for the analysis, we first built a linear model using the gene expression levels for each gene in the function of the covariates, then we selected a subset of potentially relevant genes based on their residual variance (the variability of their expression level) since genes for which there is no residual variance once the covariate effects are removed are not going to reveal any detectable effect that could be partly explained by changes at the genetic level. This approach has similarly been used in a QTL research paper for genomic methylation QTL (meQTL, [11]).
This step of identification of eQTLs confirmed a large part of previously identified couples eSNP-eGenes discovered in previous studies and identified new ones, presenting a strategy that can be applied with a limited sample size (<100).

In a second step, we aimed to identify genes that are modified during complicated pregnancy by the disease specifically and which present a component of genetic regulation. Thus, we removed the disease status from the group of covariates and submitted the data for analysis again in order to identify a novel list of eQTLs. This approach led to the discovery of 16 cis-eQTLs that are potentially associated with placental diseases. We surmise that these genetic factors identified may have been overlooked in previous analyses since they could have largely been confounded with the covariates. In particular, in our approach, we managed to partially account for critical confounding covariates such as placental age, which is often difficult to dissociate from placental disease. Interestingly, both approaches identified genes where we could observe statistical associations between disease status and genotype at the locus.

## 2.2. Human Placental Samples

This study included the use of human placental samples from three different cohorts for a total of 66 samples. For all the cohorts, the participant gave their informed consent. For the St George's cohort, the administrative ethical references were 19/LO/0974, approved by the Brent Research Ethic Committee (London) on 28 June 2019. For the Angers cohort, the study was approved by the Ethics Committee of Angers. The cohort was registered at the French CNIL (Commission Nationale de l'Informatique et des Libertés no. pWP03752UL, ethics committee for the collection of clinical data from patient records). The study was validated by the French CPP (Comité de Protection des Personnes) and registered to the French Ministry of Research under number DC-2009-907. Finally, for the Co-chin hospital cohort, the ethical administration was given by the CPP 'Ile de France XI' under the reference number 11018, 3 March 2011.

Angers University Hospital and Cochin Hospital cohorts have been described in [12]. They included 8 control (CTRL, placentas obtained from uncomplicated pregnancies) and 13 IUGR samples from Angers, 9 CTRL, 7 PE, and 3 PE+IUGR from Cochin. Preeclampsia was defined by the presence of hypertension (systolic pressure > 140 mm Hg, diastolic pressure > 90 mm Hg) and proteinuria (>0.3 g/day). IUGR was defined as a reduction of fetal growth during gestation accompanied by Doppler abnormalities and birth weight below the 10th percentile [13,14]. The Gynecochin cohort samples were kindly shared by Dr. Louis Marcellin, and included 20 CTRL, 13 from natural delivery, and 7 from C-section. For placental sample collection relative to Angers, Cochin, and Gynecochin cohorts, sections of 1 cm3 of placental villi were dissected from four different cotyledons between the basal and chorionic plates, as previously described in [14]. St. George's Hospital samples, 6 PE placentas from C-sections, were kindly provided by Prof. Basky Thilaganathan and MD Veronica Giorgione in the context of the iPlacenta consortium. PE was defined by the presence of hypertension according to guidelines by The International Society for the Study of Hypertension in Pregnancy [15]. Placental samples were collected from the middle region of the placenta containing only villous tissue, just above the cord insertion point, then washed in PBS.

## 2.3. RNA and DNA Extraction

Human placental tissues were obtained from three different cotyledons of the maternal side of the placenta, washed extensively in sterile PBS, and snap frozen or in RNAlater (Invitrogen, Waltham, MA, USA), then powdered using pestle and mortar on dry ice or the hammering method in liquid nitrogen. The powdered tissues were then processed for DNA and RNA extraction. RNA was extracted with TRIzol (Invitrogen), according to manufacturer's instructions; UltraPure™ Phenol:Chloroform:Isoamyl Alcohol (25:24:1, v/v) (Invitrogen) was used in place of Chloroform. Extracted RNA was resuspended in RNAse and DNase-free water, and integrity was assessed with Bioanalyzer Agilent 2100 nano kit.
For DNA extraction, powdered tissue was resuspended in Sample Lysis Buffer (10 mM Tris Hcl pH8, 10 mM EDTA, 50 mM NaCl, 0.5% SDS, Proteinase K (Invitrogen #25530-049) used at 1:3000) and incubated overnight in shaking water bath at 58 °C. Addition of absolute ethanol allowed precipitation of the DNA in filamentous, visible form, which was then transferred to new 1.5 mL tubes for additional washes in ethanol 70%. The pellet was left to air-dry, and DNA was resuspended in Tris-EDTA.

## 2.4. Transcriptomic Dataset

100 ng of RNA from the human placental samples were analyzed by ClariomD arrays (Applied Biosystems™, Affymetrix, Thermo Fisher Scientific, Montigny-le-Bretonneux, France) as described in [12]. This array measures gene expression both at the gene level as well as the exon level, providing also splicing isoform-specific data. The 66 samples were processed in three different batches, generating three sets of raw data files (.CEL). *Merging datasets.* To reduce batch effects in fluorescence signal due to experimental variability, the raw data (.CEL) files were processed together using the Transcriptome Analysis Console (TAC) software (Thermo Fisher Scientific), performing the default Robust Multi-Array normalization (RMA). This algorithm performs background adjustment, followed by quantile normalization and summarization [16]. The ClariomD array measures expression levels for a total of 134,748 probes, corresponding to coding and non-coding genes, 18,858 and 66,845, respectively, as well as predicted genes, pseudogenes, and small RNAs. Probe coordinates refer to the GRCh38 Genome Build. For downstream analyses, only genes with known GeneID and Description in TAC have been kept, removing transcripts identified by the database AceView, as well as transcripts on the Y and mitochondrial chromosomes; for a final total of 46,624 probes. Probes have been further filtered by mean fluorescent value (LOG2) across all samples ≥ 4.5, reducing the total number of probes to 33,988 considered to be "expressed". *Transcriptomic covariates.* Classical clinical and experimental variables in placental studies were used for the whole sample and included a total of 9 variables (Batch, Cohort, Group, Delivery, Maternal Age, Ethnicity, Gestational Age, Sex, and Parity). In particular, Batch refers to the experimental batch of transcriptomic data acquisition, Cohort to the different cohort of origins of the placental samples relative to the Hospital in which they had been collected, Group to the disease status, Delivery to the delivery mode, either C-section or Natural delivery, and finally Sex to the placental sex. For these 9 variables, missing values were replaced with the mean of each disease group, according to the disease group of the sample with missing value [17]. Summary statistics by disease group are listed in Table 1, the complete dataset is available in Supplementary Table S2.

<table><caption>**Table 1.** Summary statistics of human placental samples sorted by disease.</caption><thead><tr><th>Disease Group</th><th>CONTROLS</th><th>PE</th><th>PE + IUGR</th><th>IUGR</th></tr></thead><tbody><tr><td>Delivery mode (Caesarian/Natural)</td><td>23/12 (61.4%)</td><td>9/0 (15.8%)</td><td>3/0 (5.3%)</td><td>10/0 (17.5%)</td></tr><tr><td>Maternal age (years)</td><td>34.0 ± 3.9</td><td>34.2 ± 6.0</td><td>35.3 ± 2.4</td><td>32.1 ± 6.6</td></tr><tr><td>Ethnicity (Afr/Eur)</td><td>8/27</td><td>5/4</td><td>1/2</td><td>1/9</td></tr><tr><td>Gestational age (years)</td><td>39.2 ± 1.2</td><td>34.9 ± 2.6</td><td>30.0 ± 2.5</td><td>31.0 ± 2.8</td></tr><tr><td>Sex (M/F)</td><td>17/18</td><td>4/5</td><td>2/1</td><td>3/7</td></tr><tr><td>Parity</td><td>1.9 ± 1.5</td><td>1.2 ± 1.1</td><td>2.0 ± 1.4</td><td>1.2 ± 0.4</td></tr></tbody></table>

Principal Component Analysis (PCA) was performed in R (version 4.1.0, 2021-05-18) with the R package PCATools, on the 57 samples, for the transcriptomic dataset of 46,624 probes. The eigencorplot function within the package was used to calculate correlation coefficients between the first 10 principal components (PCs) and the known covariates; it performs a Pearson correlation, followed by F-statistic [18]. The R package SVA was used to identify potential additional sources of variation in the dataset, i.e., surrogate variables,
by building the mod linear distribution with “Group” as the independent variable and
setting the mod0, to the intercept [10]. For more details, see [12]. A correlation matrix was
calculated in R, using as input the 9 clinical and experimental variables, the first 10 PCs
of gene expression PCA, and 4 identified surrogate variables to identify the final set of
variables correlated with global transcriptomic changes. A correlation coefficient cut-off
threshold of |0.9| was used to remove colinear variables [19]. The final set of transcrip-
tomic covariates includes Batch, Group, Delivery, Maternal Age, Ethnicity, Gestational Age,
Sex, Parity, and Transcriptome PC1 to PC5.

## 2.5. Genotype Dataset

For the 66 samples, 200 ng of genomic DNA were genotyped using the Infinium
OmniExpress (illumina) BeadChip, which interrogates 713,407 SNPs. The raw data files
were analyzed with GenomeStudio2.0 software to retrieve the genotype dataset using
the A/B allele Illumina notation. SNP coordinates used refer to the GRCh38 Genome
Build. **Quality control.** Quality control (QC) of samples and SNPs was conducted in
PLINK1.9 [20,21]. Samples were subjected to quality control, removing samples with a rate
of missing genotype ≥ 2%, as well as samples that presented an excess heterozygosity as an
index of contamination F < −0.05, calculated on the pruned dataset (pairwise correlations,
window size 50 SNPs, sliding window of 5 positions, SNPs with r2 ≥ 0.2 removed).
For SNP QC, SNPs with a rate of missing genotypes across samples≥ 1%, as well as
SNPs diverging from Hardy-Weinberg equilibrium (p-value < 10⁻⁶), have been removed.
In total, 57 samples and 665,191 variants passed quality control. **Genotype covariates.**
PCA was performed in PLINK1.9 on the clean dataset (57 samples, 665,191 variants),
pruned as described above, retaining the first 10 principal components (GenotypePC1-
PC10) as covariates for downstream eQTL analysis to summarise population stratification.
**Ancestry Estimation.** Ancestry estimation has been performed to infer the ethnicity of
placental samples from Cochin and Anger cohorts. Briefly, the genotypes from the 1K
Genomes Phase III release have been downloaded from The International Genome Sample
Resource, mapped to the GRCh37 genome build, and have been used for imputation of
allele frequencies across populations [22,23]. The set of SNPs in common between the two
datasets has been established, and the two datasets merged using PLINK1.9. PCA was
performed on the merged dataset, pruned as described above, including only autosomal
chromosomes, for a total of 185,249 variants. Sample clustering along PC1 and PC2 have
been used to extrapolate the unknown ethnicities of the placental samples based on their
relative distance with the 1K Genomes Phase III samples.

## 2.6. eQTL Analysis Workflow

The R package MatrixEQTL has been used to perform the eQTL analyses in R [24], by
applying a linear multivariate model to test the contribution of each SNP genotype to gene
expression levels. Two main approaches were pursued. The ALL COVARIATES dataset
included the full set of 23 covariates (Batch, Group, Delivery Placental Sex, Gestational
Age, Maternal Age, Parity, Ethnicity, the first 5 PCs from gene expression PCA, the first
10 PCs from genotype principal component analysis that summarises population strati-
fication), allowing to identify the eSNPs correlating with changes of gene expression of
the eGenes while minimizing spurious associations between SNPs and Genes [25]. The
MINUS DISEASE dataset included all covariates but "Group" and Gene Expression PC5
given its strong correlation with the disease phenotype described by the "Group" variable.
With this approach, the contribution of the disease variable on gene expression changes
was not corrected for, allowing to investigate a set of genes for which the effects of disease
on expression are important and would have otherwise excluded to the residual variance
threshold, to identify SNPs that correlate with gene expression changes and could harbor
interaction effects with the disease phenotype.

For both approaches, the same workflow was followed. Each gene expression distribu-
tion was first normalized by a built-in function with rank ("average" method) followed by
quantile normalization to reduce the impact of outliers. The normalized gene expression dataset was then modeled by multivariate linear regression, expressing gene expression for each gene in the dataset as a function of either the ALL COVARIATES or MINUS DISEASE set of covariates. The residual variance of the model can be thought of as an indication of potential genetic influence on gene expression levels [26,27]. Using incremental thresholds of residual variance from 0.5 to 0.95 with 0.05 increments, 10 gene sets termed "RV0.05" to "RV0.95" were defined and used as input gene expression datasets in eQTL analyses.

## 2.7. *Cis-QTL Analyses of Gene Expression Subsets*

Gene subsets RV0.50 to RV0.95 and 417,114 SNPs (MAF ≥ 0.15) were analyzed for correlation with the LINEAR model parameter in MatrixEQTL [24]. We evaluated only the correlation of eGenes and eSNPs in *cis* (*cis*-eQTLs) given the limited power of our analysis and the small number of trans-QTLs previously identified in the placenta by a larger designed study [9]. To further reduce the total number of tests, the *cis*-distance has been set to 100 kb rather than the default 1 Mb [8]. The gene expression dataset has been normalized with MatrixEQTL prior to performing eQTL analysis, as described above. For the ALL COVARIATES Dataset analysis, 23 covariates (Batch, Group, Delivery, Maternal Age, Ethnicity, Gestational Age, Sex, Parity, first five PCs of gene expression PCA, first 10 PCs of genotype PCA) have been accounted for by the software when performing multi-variate linear regression for each gene-SNP pair, while for the MINUS DISEASE dataset, 21 covariates have been included: Batch, Delivery, Maternal Age, Ethnicity, Gestational Age, Sex, Parity, first four PCs of gene expression PCA, first 10 PCs of genotype PCA. To correct for multiple testing, the in-built function in MatrixEQTL has been used that estimates FDR-adjusted *p*-values (*q*-values) for each gene-SNP pair with the Benjamini-Hochberg procedure. In this study, only SNPs with FDR ≤ 0.05 have been considered as statistically significant. The full lists of statistically significant eSNP-eGene pairs for both datasets are presented in Supplementary Tables S3 and S4.

## 2.8. *Calculating Enrichment of Significant cis-QTLs for Each Subset*

To choose the best *cis*-eQTL analytical design, defined by the gene expression subset used, the enrichment of significant *cis*-QTLs has been calculated. We define enrichment as the number of *cis*-eQTL per eGene. *Cis*-eQTL enrichment has been considered together with the total number of eGenes and eSNPs identified as a result of statistically significant *cis*-eQTLs (FDR ≤ 0.05) in order to find the optimal cut-off, where we see an increase of enrichment, without losing too many eGenes. For both ALL COVARIATES and MINUS DISEASE datasets, the input gene expression dataset RV0.85 was selected as optimal for the final *cis*-QTL analysis workflow and statistically significant *cis*QTLs were further characterized and are presented here. The list of input genes with relative individual expression levels as well as residual variances, are presented in Supplementary Tables S5 and S6, respectively.

## 2.9. *Calculating Overlap with Previous Studies*

The statistical significance of the overlap with previous studies was calculated with Fisher's exact test, with normal approximation, with a statistical significance threshold set at *p*-value ≤ 0.05 [28]. As reference total number of genes, we used N = 30,818; this number of placental expressed genes, with ≥0.1 RPKM (Read Per Kilobase of exon model per Million mapped reads), as defined in the latest high-quality study on placental transcriptome by Gong and coworkers [29]. The group carried out an RNA-sequencing analysis on 302 human placental samples, including messenger RNAs, long non-coding RNAs, as well as small and circular non-coding RNAs. This number is consistent with our study, where we observe 33,988 genes with mean fluorescent value (LOG₂) ≥ 4.5. The total number of genes in the human genome, even though still lacking consensus, is reported to be higher than 35000, including coding and non-coding genes, as discussed in [30].
## 2.10. Calculating Statistical Significance of Interaction between Best-eSNP and Disease on eGene Gene Expression

To check whether an interaction exists between the genotype at the best-eSNP and the disease status on gene expression levels of the eGene we performed a linear regression in R and considered the linear model to be significant when the model *p*-value ≤ 0.05, similarly the effects of each variable on the model were considered significant when the coefficient *p*-value in the model was ≤ 0.05. The general formula of the linear regression model was: eGene ~ eSNP + Group + eSNP*Group. The linear regression was performed on the residuals of the gene expression dataset after multilinear regression with the MINUS DISEASE set of covariates to correct for confounding factors (the data had also initially been normalized by ranking and quantile normalization as described above).

# 3. Results

## 3.1. Transcriptome Identification of Confounding Variables

The transcriptome dataset is composed of 57 RNA samples that passed the quality controls from either the Cochin Hospital (Technologic Facility Gernom'IC), Angers University Hospital, or St George's Hospital, University of London (Table 1). The transcriptome data were obtained using the Clariom D microarray (Affymetrix), which allows to analysis of the level of mRNAs at the exon level [12].

From the transcriptomic datasets, only genes with known GeneID and Description in the TAC (Transcriptome Analysis Console, Affymetrix) were kept, removing transcripts identified by the database AceView, as well as transcripts on the Y and mitochondrial chromosomes (n = 616 and 18, respectively), for a final total of 46,624 probes. Probes were further filtered by mean fluorescent value (LOG2) across all samples ≥ 4.5, eventually reducing the total number of probes to 33,988.

A new analysis strategy of our work was the removal of the covariates affecting gene expression; while some are interesting biological factors that could deserve interest, such as maternal age for instance, we decided to consider them all as confounders to identify solely eQTL that influence placental gene expression without obvious medical consequences. There were nine variables available for each sample that was included ab initio in the analysis (Batch, Cohort (St Georges, Cochin, or Angers), Disease Group, Mode of Delivery, Maternal Age, Ethnicity, Gestational Age at Delivery, Sex, Parity). Using PCA in R (PCATools) enabled us to define principal components (PCA axes). The 10 first PCA axes captured ~61% of the variability, while the two first axes, PC1 and PC2 only, concentrated ~34.5% of the variability. The correlation of the PC axes with the clinical variables is represented in Figure 2A. The weights of the clinical covariates could be estimated at 7.4%, 5%, 10.7%, 5.7%, 5.8%, 14.4%, 11.1%, 23.4%, and 23.5% for Parity, Sex, Gestational age, Ethnicity, Maternal Age, Delivery Mode, Disease Group, Cohort, and Batch, respectively. The total is above 100% since there are complex correlations between the different items; for instance, disease status is correlated with gestational age. Some axes were strongly correlated with the Group (i.e., the disease status), especially PC1 and PC5, the latter one being also intimately correlated with the gestational age.

In addition, The R package SVA was used to identify possible surrogate variables (additional unknown sources of variation in the dataset). We collated all variables in a correlation matrix that was calculated in R, using as input the nine clinical and experimental variables, the first 10 PCs of gene expression PCA, and four identified surrogate variables to identify the final set of variables correlated with global transcriptomic changes. A correlation coefficient cut-off threshold of |0.9| was used to remove colinear variables [19]; see Figure 2B. For the final choice of covariates included in the eQTL analysis, we kept the five first PCs, "Batch", "Group", "Delivery", "Maternal Age", "Ethnicity", "Gestational Age" and "Sex". The covariate "Cohort" was removed since it was colinear with "Batch".
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Definition of covariates and surrogate variables potentially affecting the gene expression dataset. (A) Eigencorplot expresses the results of Pearson correlation between the clinical and technical variables and the first 10 principal components obtained from the PCA performed on the gene expression dataset. A negative correlation is expressed in blue, and a positive correlation is expressed in red. Within each pair the correlation coefficient is displayed with statistical significance. (* = *p*-value ≤ 0.05, ** = *p*-value ≤ 0.01, *** = *p*-value ≤ 0.001). (B) Correlation matrix between the 9 clinical and experimental variables, first 10 PCs of gene expression PCA, and 4 identified surrogate variables, having effects on global transcriptomic changes. Only statistically significant correlations (*p*-value ≤ 0.05) are displayed as dots of increasing size as a measure of the correlation coefficient. Positive correlations are displayed in blue, and negative correlations are displayed in red. In yellow boxes are the variables that were kept as covariates for subtracting their effect (see text).</figcaption></figure>

## 3.2. Genotyping and Population Stratification

The final genotype dataset included 665,191 variants for the 57 samples, of which we kept SNPs with a Minor Allele Frequency > 0.15 (417,114 were kept). The legislation in France forbids the collection of ethnic data. Nevertheless, with the genotype information, we were able to classify our samples using PCA alone and with reference to the 1000 genome project (Supplementary Figure S1). Seven placental samples were clearly of African origin, 42 were of European origin, and eight were mixtures between African and European backgrounds. The first 10 PCs from the genotype principal component analysis were included as covariates in the eQTL analysis to account for variations due only to population stratification.

## 3.3. Optimal Feature Selection for the eQTL Analysis

We first explored the effect of varying the number of input genes for the eQTL analysis to identify the optimal gene sets in order to achieve a satisfactory balance between the number of statistically significant cis-eQTLs (FDR ≤ 0.05), the number of identified eGenes and cis-eQTL enrichment (in 100 kb proximal to the gene), defined as the number of eSNPs/eGene. For this, we used the concept of residual variance (RV) to reduce the number of tested genes.
Given our limited sample size, it was our primary goal to be able to define the set of genes to be used as input in the eQTL analysis that would give us the optimal experimental design to identify significant cis-eQTLs, reducing as much as possible the number of tests, without missing out potentially relevant genes. The eQTL workflow that we used was based on multivariate linear modeling, which expresses the variance observed in the expression levels for each gene in the function of the input covariates. Therefore, the RV is a measure of potential genetic influence on gene expression levels, amongst other variables [26,27]. The summary of the threshold optimization for analysis is presented in Supplementary Table S1. While the number of significant eGenes decreased from 153 to 22, when the RV raised from 0.5 to 0.95, the number of significant SNP per gene increased from 4.00 to 8.91.

We selected an RV threshold of 0.85 for further analysis as a good compromise between having significant but relevant eQTL, albeit other thresholds could certainly be chosen. At this threshold, at least six SNP per gene were significant, meaning that the association with gene expression is not merely due to linkage disequilibrium. The analysis at this threshold rested upon 3201 genes, 417114 SNPs, a multivariate model with 23 covariates, and a window of ±100 kb around each gene, following the experimental choice of Kikas et al. [8], which justified this window by the fact that relevant regulatory SNPs are mostly in the vicinity (~100 kb) of the gene. These parameters resulted in the discovery of 279 statistically significant placental SNPs (p < 1.24 × 10⁻⁴, FDR < 0.05) proximal to 43 eGenes (Figure 3, Table 2). The full list of statistically significant eGene-eSNP pairs is provided in Supplementary Table S3.

In Table 2, the novel genes discovered are presented in yellow, while the genes in bold correspond to the most stringent threshold and to a Bernoulli genome-wide correction which does not consider the possible LD relating to successive SNPs. When this is considered, the Bernoulli correction corresponds to a limited number of independent tests, leading to a threshold of ~10⁻⁴ instead of 10⁻⁸ (blue and red line, respectively, in Figure 3A). We report 22 novel placental eQTL, while 51% of the placental eGenes identified here have been previously described as placental cis-eQTLs in the three placental studies available to date (n = seventeen for the Peng et al. dataset, n = eight for Delahaye et al., and n = five for Kikas et al.) [7–9]. Fisher’s exact tests of the overlaps between datasets were all statistically significant with *p*-values < 0.05, respectively *p* < 5.83 × 10⁻⁷, *p* < 1.88 × 10⁻⁶, and *p* < 2.75 × 10⁻⁸, validating our experimental approach.

Interestingly, out of the 279 cis-QTL found, 37 corresponded to SNPs located between two genes modified at the expression level. On chromosome 6, 24 eSNPs were found to be shared between TOB2P1 and ZSCAN9 out of their 41 and 42 cis-eQTLs, respectively (Figure 3B). On chromosome 11, two different regions of interest were identified, between LOC646029 and AQP11, that shared 10 eSNPs, out of 11 and 12 cis eQTLs, respectively. Finally, three eSNPs associated with CARD17 expression deregulations correlated also with upstream transcript CASP1P2, which had a total of 11 cis-eQTLs (Figure 3B). These observations strongly suggest that the SNPs influence gene expression in the chromosome region, either by modifying an enhancer acting on two genes close by or by being transcribed in a non-coding RNA having a cis-regulatory effect on transcription.
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Overview of the cisQTL results for the ALL COVARIATES dataset. (A) Manhattan plot showing each tested SNP as a dot in its genomic position on the x-axis, stratified by chromosomes. The y-axis expresses the -LOG10 of the *p*-value of the interaction with the eGene. Suggestive line in red shows the $10 \times 10^{-8}$ threshold, while in blue, the threshold corresponding to a *p*-value of 0.05 adjusted for the number of tests (False Discovery Rate) = $1.24 \times 10^{-4}$. (B) Plots showing the localization of the tested SNPs within the genomic regions on chromosome 6, containing ZSCAN9 and TOB2P1, and chromosome 11, where AQP11, LOC646029, CASP1P2, and CARD17 are located. The x-axis expresses the genomic coordinates, while in the y-axis, the -LOG10 of the *p*-value, the horizontal line sets the level for FDR ≤ 0.05. The vertical lines represent the start and end sites for each of the genes. eSNPs that were found to be in common between eGenes in the same region are highlighted in yellow.</figcaption></figure>
<table><caption>**Table 2. List of eGenes, with the most significant cis-eQTL for the ALL COVARIATES dataset.** The table lists the total number of eGenes correlating with the identified significant cis-eQTLs; for each eGene the total number of eSNPs are listed, and the details of the Best-eSNP reported. The Best-eSNP, in this case, is defined as the most significant cis-eQTL associated with each eGene (smallest recorded *q*-value). Beta corresponds to the slope coefficient for the multivariate linear model and is an estimation of the eSNP effect size on gene expression variance. Novel placental eGenes are highlighted in yellow. If the eGene had been previously identified as such, the studies are specified in the “Previous study” column.</caption><thead><tr><th>Gene Symbol</th><th>Description</th><th>eSNPs</th><th>Best-eSNP</th><th>p-Value</th><th>q-Value</th><th>Beta</th><th>Chr</th><th>Pos.</th><th>Previous Study</th></tr></thead><tbody><tr><td>ZSCAN9</td><td>zinc finger and SCAN domain containing 9</td><td>42</td><td>rs1150707</td><td>6.46 × 10<sup>-16</sup></td><td>1.43 × 10<sup>-11</sup></td><td>18.80</td><td>6</td><td>28229827</td><td>[8]</td></tr><tr><td>PSG7</td><td>pregnancy specific beta-1-glycoprotein 7 (gene/pseudogene)</td><td>3</td><td>rs7248225</td><td>3.20 × 10<sup>-12</sup></td><td>1.68 × 10<sup>-8</sup></td><td>27.40</td><td>19</td><td>42918847</td><td>[7-9]</td></tr><tr><td>TOB2P1</td><td>transducer of ERBB2, 2 pseudogene 1</td><td>41</td><td>rs13408</td><td>4.79 × 10<sup>-12</sup></td><td>2.21 × 10<sup>-8</sup></td><td>-17.74</td><td>6</td><td>28244970</td><td>[9]</td></tr><tr><td>ERAP2</td><td>endoplasmic reticulum aminopeptidase 2</td><td>13</td><td>rs2549778</td><td>8.01 × 10<sup>-12</sup></td><td>3.28 × 10<sup>-8</sup></td><td>19.34</td><td>5</td><td>96868551</td><td>[7-9]</td></tr><tr><td>AQP11</td><td>aquaporin 11</td><td>12</td><td>rs10793257</td><td>4.80 × 10<sup>-11</sup></td><td>1.29 × 10<sup>-7</sup></td><td>20.84</td><td>11</td><td>77598488</td><td>[7-9]</td></tr><tr><td>LOC646029</td><td>uncharacterized LOC646029</td><td>11</td><td>rs10793257</td><td>2.48 × 10<sup>-10</sup></td><td>6.09 × 10<sup>-7</sup></td><td>20.61</td><td>11</td><td>77598488</td><td></td></tr><tr><td>RPS26</td><td>ribosomal protein S26</td><td>9</td><td>rs11171739</td><td>8.12 × 10<sup>-10</sup></td><td>1.95 × 10<sup>-6</sup></td><td>18.98</td><td>12</td><td>56076841</td><td>[7]</td></tr><tr><td>RP1-97J1.2</td><td>putative novel transcript</td><td>7</td><td>rs9372316</td><td>1.14 × 10<sup>-9</sup></td><td>2.62 × 10<sup>-6</sup></td><td>21.02</td><td>6</td><td>112000000</td><td></td></tr><tr><td>MLLT10</td><td>myeloid/lymphoid or mixed-lineage leukemia</td><td>15</td><td>rs10828248</td><td>3.35 × 10<sup>-9</sup></td><td>6.73 × 10<sup>-6</sup></td><td>-16.45</td><td>10</td><td>21535690</td><td></td></tr><tr><td>IL36RN</td><td>interleukin 36 receptor antagonist</td><td>25</td><td>rs6761276</td><td>7.09 × 10<sup>-9</sup></td><td>1.31 × 10<sup>-5</sup></td><td>-17.30</td><td>2</td><td>113000000</td><td>[8,9]</td></tr><tr><td>LYNX1</td><td>Ly6/neurotoxin 1</td><td>9</td><td>rs10956986</td><td>1.13 × 10<sup>-7</sup></td><td>1.47 × 10<sup>-4</sup></td><td>16.92</td><td>8</td><td>143000000</td><td>[9]</td></tr><tr><td>CASP1P2</td><td>caspase 1 pseudogene 2</td><td>11</td><td>rs1023954</td><td>1.18 × 10<sup>-7</sup></td><td>1.50 × 10<sup>-4</sup></td><td>19.23</td><td>11</td><td>105000000</td><td></td></tr><tr><td>MIR4804</td><td>microRNA 4804</td><td>5</td><td>rs2253215</td><td>1.26 × 10<sup>-7</sup></td><td>1.56 × 10<sup>-4</sup></td><td>23.89</td><td>5</td><td>72952041</td><td></td></tr><tr><td>ZFYVE19</td><td>zinc finger, FYVE domain containing 19</td><td>6</td><td>rs10152371</td><td>2.43 × 10<sup>-7</sup></td><td>2.63 × 10<sup>-4</sup></td><td>16.13</td><td>15</td><td>40811095</td><td>[7,9]</td></tr><tr><td>HJURP</td><td>Holliday junction recognition protein</td><td>4</td><td>rs2361506</td><td>4.05 × 10<sup>-7</sup></td><td>4.11 × 10<sup>-4</sup></td><td>18.92</td><td>2</td><td>234000000</td><td></td></tr><tr><td>ANAPC4</td><td>anaphase promoting complex subunit 4</td><td>10</td><td>rs1993602</td><td>5.38 × 10<sup>-7</sup></td><td>5.41 × 10<sup>-7</sup></td><td>-15.90</td><td>4</td><td>25413484</td><td>[9]</td></tr><tr><td>OR2F1</td><td>olfactory receptor, family 2, subfamily F, member 1 (gene/pseudogene)</td><td>5</td><td>rs7798409</td><td>1.08 × 10<sup>-6</sup></td><td>1.00 × 10<sup>-3</sup></td><td>-20.39</td><td>7</td><td>144000000</td><td></td></tr><tr><td>MIR3927</td><td>microRNA 3927</td><td>1</td><td>rs7046565</td><td>1.72 × 10<sup>-6</sup></td><td>1.56 × 10<sup>-3</sup></td><td>16.92</td><td>9</td><td>110000000</td><td></td></tr><tr><td>GBP3</td><td>guanylate binding protein 3</td><td>5</td><td>rs12121223</td><td>2.50 × 10<sup>-6</sup></td><td>2.14 × 10<sup>-3</sup></td><td>-18.83</td><td>1</td><td>89015900</td><td>[9]</td></tr><tr><td>FUT10</td><td>fucosyltransferase 10 (alpha (1,3) fucosyltransferase)</td><td>5</td><td>rs7018447</td><td>3.31 × 10<sup>-6</sup></td><td>2.67 × 10<sup>-3</sup></td><td>16.61</td><td>8</td><td>33467537</td><td>[9]</td></tr><tr><td>NSUN7</td><td>NOP2/Sun domain family, member 7</td><td>3</td><td>rs2437317</td><td>6.06 × 10<sup>-6</sup></td><td>4.29 × 10<sup>-3</sup></td><td>25.15</td><td>4</td><td>40789990</td><td></td></tr><tr><td>CBLB</td><td>Cbl proto-oncogene B, E3 ubiquitin protein ligase</td><td>1</td><td>rs1503920</td><td>6.32 × 10<sup>-6</sup></td><td>4.45 × 10<sup>-3</sup></td><td>15.58</td><td>3</td><td>106000000</td><td>[8,9]</td></tr><tr><td>AC104135.2</td><td>novel transcript</td><td>3</td><td>rs7573356</td><td>9.72 × 10<sup>-6</sup></td><td>6.43 × 10<sup>-3</sup></td><td>-20.13</td><td>2</td><td>74941537</td><td></td></tr></tbody></table>
<table><caption>**Table 2. Cont.**</caption><thead><tr><th>Gene Symbol</th><th>Description</th><th>eSNPs</th><th>Best-eSNP</th><th>p-Value</th><th>q-Value</th><th>Beta</th><th>Chr</th><th>Pos.</th><th>Previous Study</th></tr></thead><tbody><tr><td>GPR132</td><td>G protein-coupled receptor 132</td><td>2</td><td>rs7157567</td><td>1.05 × 10<sup>-5</sup></td><td>6.70 × 10<sup>-3</sup></td><td>18.27</td><td>14</td><td>105000000</td><td></td></tr><tr><td>GTSF1</td><td>gametocyte specific factor 1</td><td>1</td><td>rs11170917</td><td>1.31 × 10<sup>-5</sup></td><td>7.98 × 10<sup>-3</sup></td><td>24.42</td><td>12</td><td>54472204</td><td>[9]</td></tr><tr><td>SLC13A5</td><td>solute carrier family 13 (sodium-dependent citrate transporter), member 5</td><td>1</td><td>rs9889374</td><td>1.46 × 10<sup>-5</sup></td><td>8.69 × 10<sup>-3</sup></td><td>17.28</td><td>17</td><td>6662428</td><td></td></tr><tr><td>CASQ1</td><td>calsequestrin 1 (fast-twitch, skeletal muscle)</td><td>1</td><td>rs6693877</td><td>1.67 × 10<sup>-5</sup></td><td>9.80 × 10<sup>-3</sup></td><td>-14.62</td><td>1</td><td>160000000</td><td>[9]</td></tr><tr><td>AC009542.2</td><td>novel transcript, antisense to WDR91</td><td>3</td><td>rs10954493</td><td>1.87 × 10<sup>-5</sup></td><td>1.08 × 10<sup>-2</sup></td><td>-18.52</td><td>7</td><td>135000000</td><td></td></tr><tr><td>KANSL1-AS1</td><td>KANSL1 antisense RNA 1</td><td>1</td><td>rs17585974</td><td>2.62 × 10<sup>-5</sup></td><td>1.46 × 10<sup>-2</sup></td><td>24.81</td><td>17</td><td>46171833</td><td></td></tr><tr><td>MTPAP</td><td>mitochondrial poly(A) polymerase</td><td>2</td><td>rs1762598</td><td>4.17 × 10<sup>-5</sup></td><td>2.15 × 10<sup>-2</sup></td><td>-15.63</td><td>10</td><td>30352501</td><td></td></tr><tr><td>CPS1</td><td>carbamoyl-phosphate synthase 1</td><td>3</td><td>rs918233</td><td>4.27 × 10<sup>-5</sup></td><td>2.20 × 10<sup>-2</sup></td><td>-16.64</td><td>2</td><td>211000000</td><td>[9]</td></tr><tr><td>MIR4299</td><td>microRNA 4299</td><td>1</td><td>rs7126296</td><td>4.53 × 10<sup>-5</sup></td><td>2.32 × 10<sup>-2</sup></td><td>17.44</td><td>11</td><td>11556715</td><td></td></tr><tr><td>AGA</td><td>aspartylglucosaminidase</td><td>2</td><td>rs4690523</td><td>5.31 × 10<sup>-5</sup></td><td>2.63 × 10<sup>-2</sup></td><td>-15.80</td><td>4</td><td>177000000</td><td>[9]</td></tr><tr><td>SOHLH2</td><td>spermatogenesis and oogenesis specific basic helix-loop-helix 2</td><td>4</td><td>rs9575600</td><td>5.41 × 10<sup>-5</sup></td><td>2.64 × 10<sup>-2</sup></td><td>16.53</td><td>13</td><td>36223006</td><td>[7]</td></tr><tr><td>C8orf46</td><td>chromosome 8 open reading frame 46</td><td>1</td><td>rs12541098</td><td>6.17 × 10<sup>-5</sup></td><td>2.90 × 10<sup>-2</sup></td><td>-15.38</td><td>8</td><td>66424147</td><td>[7-9]</td></tr><tr><td>CTC-482H14.5</td><td>novel transcript, antisense to PTPRS</td><td>1</td><td>rs2251127</td><td>6.29 × 10<sup>-5</sup></td><td>2.90 × 10<sup>-2</sup></td><td>16.34</td><td>19</td><td>5138218</td><td></td></tr><tr><td>NDUFS5</td><td>NADH dehydrogenase (ubiquinone) Fe-S protein 5, 15kDa (NADH-coenzyme Q reductase)</td><td>2</td><td>rs10888650</td><td>6.30 × 10<sup>-5</sup></td><td>2.90 × 10<sup>-2</sup></td><td>-13.27</td><td>1</td><td>39041489</td><td>[8,9]</td></tr><tr><td>VN1R84P</td><td>vomeronasal 1 receptor 84 pseudogene</td><td>1</td><td>rs2015481</td><td>8.31 × 10<sup>-5</sup></td><td>3.61 × 10<sup>-2</sup></td><td>13.59</td><td>19</td><td>21676192</td><td></td></tr><tr><td>LOC101928775</td><td>uncharacterized LOC101928775</td><td>1</td><td>rs10982832</td><td>8.33 × 10<sup>-5</sup></td><td>3.61 × 10<sup>-2</sup></td><td>16.77</td><td>9</td><td>116000000</td><td></td></tr><tr><td>FGF19</td><td>fibroblast growth factor 19</td><td>1</td><td>rs7105655</td><td>1.01 × 10<sup>-4</sup></td><td>4.24 × 10<sup>-2</sup></td><td>13.69</td><td>11</td><td>69738836</td><td></td></tr><tr><td>CARD17</td><td>caspase recruitment domain family, member 17</td><td>3</td><td>rs1623342</td><td>1.09 × 10<sup>-4</sup></td><td>4.42 × 10<sup>-2</sup></td><td>15.14</td><td>11</td><td>105000000</td><td>[7,9]</td></tr><tr><td>MIR4527</td><td>microRNA 4527</td><td>1</td><td>rs982265</td><td>1.17 × 10<sup>-4</sup></td><td>4.67 × 10<sup>-2</sup></td><td>-16.51</td><td>18</td><td>47444531</td><td></td></tr><tr><td>Y RNA</td><td>Y RNA</td><td>1</td><td>rs2248978</td><td>1.18 × 10<sup>-4</sup></td><td>4.69 × 10<sup>-2</sup></td><td>13.64</td><td>12</td><td>105000000</td><td></td></tr></tbody></table>

Highlighted in light yellow are the genes that were not found in the previous studies.
### 3.4. Identifying eQTLs Involved in the Disease by a Subtraction Strategy

The same cis-eQTL analysis was performed de novo, suppressing the parameter 'Disease' from the confounding factors, as well as PC5 (Figure 2), which were very highly correlated. This allowed us to include in the analysis genes for which the disease has a strong influence on changes in gene expression levels, reducing the overall residual variance after modeling (the number of genes increased from 3201 to 3279, Figure 4). However, these genes are also believed to be the most interesting in terms of understanding pregnancy pathologies. For this reason, we aimed to investigate the contribution of genotype changes on these genes in the context of disease to be able to also detect potential genotype-disease interactions that would suggest the presence of conditional eQTLs [31] Interestingly, in this case, the number of identified eGenes was always higher than when the 'Disease' was used as a covariate. At the chosen threshold of 0.85 for the residual variance (Figure 4), the number of eGenes identified increases from 43 to 54, with the number of statistically significant (FDR ≤0.05) eSNP increasing from 279 to 307, as represented in Figure 4. The interpretation of this figure is that 16 SNPs are specific only to the disease, 38 are influenced by the covariates and the disease, while amongst the five specific to the dataset analyzed with all the covariates, four did not pass the RV threshold, MLLT10, HJURP, MIR4527, and Y_RNA. The last one, FGF19, is present in both datasets, but when the disease is not considered as a covariate, FGF19 is no longer observed as an eGene. The full list of statistically significant eGene-eSNP pairs is presented in Supplementary Table S4, while the top-eSNP is listed in Supplementary Table S7.

<figure><img src="image_7.png" /><figcaption>**Figure 4.** Comparing the experimental parameters for the ALL COVARIATES and MINUS DISEASE eQTL datasets. (TOP) Table comparing the parameters of cis-eQTL analyses performed for different Residual Variance (RV) thresholds (0.5-0.95) for ALL COVARIATES and MINUS DISEASE datasets. Depending on the RV threshold, different lists of genes were included in the cisQTL analysis, resulting in 10 output cis-QTL datasets for each approach, ALL COVARIATES in grey and MINUS DISEASE in white. (BOTTOM) Cis-eQTL enrichment graphs for each approach at the different thresholds of residual variance tested. Showing the number of significant eGenes (blue), eSNP (orange), and enrichment expressed as number of eSNP/eGene in grey.</figcaption></figure>
The 16 eGenes that were specific to the disease are presented in Figure 5 and Table 3. Several analyses of gene ontology were carried out using these 16 eGenes (disease-specific) or the 54 eGenes (placental-specific), including the 38 genes that were also presented in the first analysis using String (https://string-db.org/, accessed on 1 July 2021). All these analyses failed to identify a significant grouping of genes, suggesting that the couples eSNP-eGene that we could identify correspond to various aspects of placental biology and physiopathology and not to several specific categories or gene ontologies. Not surprisingly, however, some of these genes have been previously identified as deregulated in pregnancy pathologies (Table 3).

<figure><img src="image_8.png" /><figcaption>**Figure 5.** Overlap between cis-QTL datasets ALL COVARIATES and MINUS DISEASE. Venn diagrams showing the overlap between the eSNPs and eGenes identified as statistically significant (FDR ≤ 0.05) in ALL COVARIATES (blue) and MINUS DISEASE (yellow).</figcaption></figure>
<table><caption>**Table 3.** The 16 placental disease-specific eQTL.</caption><thead><tr><th>Gene Symbol</th><th>Description</th><th>Chr</th><th>Strand</th><th>Group</th><th>Keyword</th><th>References</th></tr></thead><tbody><tr><td>ERICH1</td><td>glutamate rich 1</td><td>chr8</td><td>-</td><td>Multiple_Complex</td><td>PE, IUGR</td><td>[32]<br>Identified as eGene by [9]</td></tr><tr><td>DNAJC15</td><td>DnaJ (Hsp40) homolog, subfamily C, member 15</td><td>chr13</td><td>+</td><td>Multiple_Complex</td><td>PE</td><td>Identified as eGene by [8,9]</td></tr><tr><td>TAS2R64P</td><td>taste receptor, type 2, member 64, pseudogene</td><td>chr12</td><td>-</td><td>Multiple_Complex</td><td></td><td></td></tr><tr><td>NSF</td><td>N-ethylmaleimide-sensitive factor</td><td>chr17</td><td>+</td><td>Multiple_Complex</td><td>Membrane fusion</td><td>[33]<br>Identified as eGene by [7,9]</td></tr><tr><td>NPM1P35</td><td>nucleophosmin 1 (nucleolar phosphoprotein B23, numatrin) pseudogene 35</td><td>chr11</td><td>+</td><td>Pseudogene</td><td>NPM1P35</td><td></td></tr><tr><td>PTTG1</td><td>pituitary tumor-transforming 1</td><td>chr5</td><td>+</td><td>Multiple_Complex</td><td>Trophoblast invasion</td><td>[34]<br>Identified as eGene by [9]</td></tr><tr><td>C8orf89</td><td>chromosome 8 open reading frame 89</td><td>chr8</td><td>-</td><td>Multiple_Complex</td><td>C8orf89</td><td></td></tr><tr><td>LINC00654</td><td>long intergenic non-protein coding RNA 654</td><td>chr20</td><td>-</td><td>NonCoding</td><td>LINC00654</td><td></td></tr><tr><td>LINC01273</td><td>long intergenic non-protein coding RNA 1273</td><td>chr20</td><td>+</td><td>NonCoding</td><td>LINC01273</td><td></td></tr><tr><td>FGF12-AS1</td><td>FGF12 antisense RNA 1</td><td>chr3</td><td>+</td><td>NonCoding</td><td>tumor suppressor</td><td>[35]</td></tr><tr><td>MYRF</td><td>myelin regulatory factor</td><td>chr11</td><td>+</td><td>Multiple_Complex</td><td>autism</td><td>[36]</td></tr><tr><td>ITPKB-IT1</td><td>ITPKB intronic transcript 1</td><td>chr1</td><td>-</td><td>NonCoding</td><td></td><td></td></tr><tr><td>CITF22-49E9.3</td><td>novel transcript</td><td>chr22</td><td>-</td><td>NonCoding</td><td></td><td></td></tr><tr><td>SLC30A8</td><td>solute carrier family 30 (zinc transporter), member 8</td><td>chr8</td><td>+</td><td>Multiple_Complex</td><td>gestational weight gain, diabetes</td><td>[37]</td></tr><tr><td>PROSER2-AS1</td><td>PROSER2 antisense RNA 1</td><td>chr10</td><td>-</td><td>NonCoding</td><td>Placental imprinted, risk for pediatric fracture,</td><td>[38]</td></tr><tr><td>PKN3</td><td>protein kinase N3</td><td>chr9</td><td>+</td><td>Multiple_Complex</td><td>endothelial cell activation, angiogenesis</td><td>[39]</td></tr></tbody></table>

-/+ stands for the reference strand of DNA on which the gene is situated in the coding orientation.

While previously identified eQTLs, such as influencing ZSCAN9 and ERAP2, have shown a strong genetic influence on their regulation, the implication of this on placental function remains elusive. In Figures 6–8, we show the levels of expression stratified according to genotype of the most significant eSNP for a subset of eGenes that presented a pronounced conditional response depending on disease group (Control vs Disease) and have been previously linked to placental function and diseases such as IL36R, FUT10, PTTG1, CBLB, and DNAJC15. Here, we plot the raw data for the expression levels; however, during analysis, the data is first normalized and then corrected for the covariates. Table 4 summarizes the results of Genotype*Disease interaction testing for each of the eGenes identified and its most significant eSNP. Interestingly, ZSCAN9, AQP11, LOC646029, FUT10, NDUFS5, and C8ORF89 showed statistically significant (p-value ≤ 0.05) interaction between the effects of disease status and the eSNP genotype on the expression of the eGene. This suggests the existence of a complex interaction between the gene and the phenotype of the disease.
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Boxplots of eGene expression in relation to the best-eSNP genotype, stratified by Group: Control vs. Disease. eGenes ZSCAN9, ERAP2, AQP11, RP1-97J1.2, IL36RN, ZFYVE19, FUT10, PTTG1, DNAJC15. The genotype of each SNP is expressed as 0, 1, 2, and the gene expression data corresponds to the raw data as mean fluorescence value. CTRL = Control; IUGR = Intra-uterine Growth Restriction; PE = Preeclampsia. Genes discussed in the text, because of their relevance to placental biology, are highlighted in green.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Boxplots of eGene expression in relation to the best-eSNP genotype, stratified by Group: Control vs. Disease. eGenes TAS2R64P, AC009542.2, CBLB, GTSF1, NDUFS5, C8orf46, LINC00654, SLC13A5, MIR4299. The genotype of each SNP is expressed as 0, 1, 2, and the gene expression data corresponds to the raw data as mean fluorescence value. CTRL = Control; IUGR = Intra-uterine Growth Restriction; PE = Preeclampsia. Genes discussed in the text, because of their relevance to placental biology, are highlighted in green.</figcaption></figure>
<table><caption>**Table 4.** Best-eSNP-eGene pairs for the MINUS DISEASE dataset showing the summary statistics from the eQTL analysis and the results of testing the significance of the interaction between genotype and disease on gene expression. The Best-eSNP, in this case, is defined as the most significant cis-eQTL associated with each eGene (smallest recorded *q*-value). Beta corresponds to the slope coefficient for the multivariate linear model and is an estimation of the eSNP effect size on gene expression variance. The table shows the test *p*-values for each coefficient of the linear model to test the effect of genotype and disease on gene expression and their interaction, as well as the model *p*-value. eGenes that presented a statistically significant interaction between eSNP genotype and disease status (*p*-value ≤ 0.05) are highlighted in green.</caption><thead><tr><th rowspan="2">eGene</th><th rowspan="2">Best-eSNP</th><th colspan="3">From MatrixEQTL Analysis MINUS DISEASE</th><th colspan="5">From Linear Regression to Test Genotype-Disease Interaction</th></tr><tr><th>p-Value</th><th>FDR</th><th>beta</th><th>Intercept p-Val</th><th>eSNP p-Val</th><th>Group p-Val</th><th>eSNP*Group p-Val</th><th>Model p-Val</th></tr></thead><tbody><tr><td>ZSCAN9</td><td>rs1150707</td><td>1.50 × 10<sup>-16</sup></td><td>3.39 × 10<sup>-12</sup></td><td>18.807</td><td>0.102647651</td><td>0.022415409</td><td>0.00787265</td><td>0.007064547</td><td>2.70 × 10<sup>-11</sup></td></tr><tr><td>ERAP2</td><td>rs2549778</td><td>1.27 × 10<sup>-12</sup></td><td>6.63 × 10<sup>-9</sup></td><td>18.894</td><td>0.237759766</td><td>0.000807833</td><td>0.029435982</td><td>0.310364841</td><td>2.02 × 10<sup>-9</sup></td></tr><tr><td>TOB2P1</td><td>rs13408</td><td>1.38 × 10<sup>-12</sup></td><td>6.63 × 10<sup>-9</sup></td><td>-17.608</td><td>0.004647281</td><td>0.001326778</td><td>0.886480707</td><td>0.50710562</td><td>1.53 × 10<sup>-8</sup></td></tr><tr><td>PSG7</td><td>rs7248225</td><td>2.47 × 10<sup>-12</sup></td><td>1.03 × 10<sup>-8</sup></td><td>26.757</td><td>0.01610927</td><td>0.000161224</td><td>0.644681149</td><td>0.935717136</td><td>1.80 × 10<sup>-8</sup></td></tr><tr><td>AQP11</td><td>rs10793257</td><td>1.10 × 10<sup>-11</sup></td><td>3.37 × 10<sup>-8</sup></td><td>20.693</td><td>0.936757201</td><td>0.244954435</td><td>0.002232625</td><td>0.009747518</td><td>3.07 × 10<sup>-8</sup></td></tr><tr><td>RP1-97J1.2</td><td>rs9372316</td><td>1.30 × 10<sup>-10</sup></td><td>3.42 × 10<sup>-7</sup></td><td>21.469</td><td>0.04677023</td><td>0.02282663</td><td>0.98324046</td><td>0.5861815</td><td>1.34 × 10<sup>-5</sup></td></tr><tr><td>LOC646029</td><td>rs10793257</td><td>1.59 × 10<sup>-10</sup></td><td>3.91 × 10<sup>-7</sup></td><td>20.487</td><td>0.62483721</td><td>0.449090131</td><td>0.000702104</td><td>0.003472929</td><td>2.96 × 10<sup>-8</sup></td></tr><tr><td>RPS26</td><td>rs11171739</td><td>3.39 × 10<sup>-10</sup></td><td>8.00 × 10<sup>-7</sup></td><td>18.460</td><td>0.05049322</td><td>0.04455226</td><td>0.32794464</td><td>0.08493583</td><td>3.95 × 10<sup>-7</sup></td></tr><tr><td>IL36RN</td><td>rs6761276</td><td>2.54 × 10<sup>-9</sup></td><td>5.23 × 10<sup>-6</sup></td><td>-17.273</td><td>0.04374326</td><td>0.00261548</td><td>0.37625035</td><td>0.56160599</td><td>3.05 × 10<sup>-7</sup></td></tr><tr><td>CASP1P2</td><td>rs1023954</td><td>3.22 × 10<sup>-8</sup></td><td>4.76 × 10<sup>-5</sup></td><td>19.003</td><td>0.08623347</td><td>0.01572294</td><td>0.91051179</td><td>0.45255771</td><td>2.25 × 10<sup>-5</sup></td></tr><tr><td>LYNX1</td><td>rs10956986</td><td>3.60 × 10<sup>-8</sup></td><td>4.86 × 10<sup>-5</sup></td><td>17.143</td><td>0.077423009</td><td>0.005369633</td><td>0.261583699</td><td>0.381727504</td><td>6.46 × 10<sup>-7</sup></td></tr><tr><td>ANAPC4</td><td>rs1993602</td><td>1.00 × 10<sup>-7</sup></td><td>1.26 × 10<sup>-4</sup></td><td>-16.357</td><td>0.62788117</td><td>0.22099995</td><td>0.06658666</td><td>0.10987541</td><td>4.00 × 10<sup>-6</sup></td></tr><tr><td>MIR4804</td><td>rs2253215</td><td>1.90 × 10<sup>-7</sup></td><td>2.15 × 10<sup>-4</sup></td><td>22.803</td><td>0.73283414</td><td>0.05996773</td><td>0.08001675</td><td>0.27408907</td><td>6.05 × 10<sup>-6</sup></td></tr><tr><td>ZFYVE19</td><td>rs10152371</td><td>2.32 × 10<sup>-7</sup></td><td>2.54 × 10<sup>-4</sup></td><td>16.012</td><td>0.29373204</td><td>0.05488449</td><td>0.18253459</td><td>0.21546712</td><td>4.07 × 10<sup>-6</sup></td></tr><tr><td>AC104135.2</td><td>rs7573356</td><td>5.00 × 10<sup>-7</sup></td><td>4.97 × 10<sup>-4</sup></td><td>-20.261</td><td>0.1301216</td><td>0.146468</td><td>0.7380213</td><td>0.2503687</td><td>3.31 × 10<sup>-4</sup></td></tr><tr><td>GBP3</td><td>rs12121223</td><td>8.00 × 10<sup>-7</sup></td><td>7.81 × 10<sup>-4</sup></td><td>-19.408</td><td>0.025397469</td><td>0.004030585</td><td>0.591077883</td><td>0.775339911</td><td>8.26 × 10<sup>-6</sup></td></tr><tr><td>OR2F1</td><td>rs7798409</td><td>9.03 × 10<sup>-7</sup></td><td>8.59 × 10<sup>-4</sup></td><td>-20.163</td><td>0.2916771</td><td>0.1675941</td><td>0.5381669</td><td>0.4819719</td><td>1.89 × 10<sup>-3</sup></td></tr><tr><td>FUT10</td><td>rs12155581</td><td>2.26 × 10<sup>-6</sup></td><td>1.93 × 10<sup>-3</sup></td><td>17.033</td><td>0.931411621</td><td>0.806002691</td><td>0.013310468</td><td>0.004690609</td><td>5.07 × 10<sup>-6</sup></td></tr></tbody></table>
<table><caption>Table 4. *Cont.*</caption><thead><tr><th rowspan="2">eGene</th><th rowspan="2">Best-eSNP</th><th colspan="3">From MatrixEQTL Analysis MINUS DISEASE</th><th colspan="5">From Linear Regression to Test Genotype-Disease Interaction</th></tr><tr><th>p-Value</th><th>FDR</th><th>beta</th><th>Intercept p-Val</th><th>eSNP p-Val</th><th>Group p-Val</th><th>eSNP*Group p-Val</th><th>Model p-Val</th></tr></thead><tbody><tr><td>MIR3927</td><td>rs7046565</td><td>2.41 × 10<sup>-6</sup></td><td>2.04 × 10<sup>-3</sup></td><td>16.732</td><td>0.020431665</td><td>0.002434153</td><td>0.405350891</td><td>0.294832332</td><td>5.07 × 10<sup>-4</sup></td></tr><tr><td>GPR132</td><td>rs7157567</td><td>3.86 × 10<sup>-6</sup></td><td>2.96 × 10<sup>-3</sup></td><td>18.710</td><td>0.06895758</td><td>0.03121023</td><td>0.83496731</td><td>0.76965232</td><td>2.64 × 10<sup>-3</sup></td></tr><tr><td>ERICH1</td><td>rs1703911</td><td>5.35 × 10<sup>-6</sup></td><td>3.90 × 10<sup>-3</sup></td><td>-13.888</td><td>0.048146347</td><td>0.003527824</td><td>0.9634337</td><td>0.940283908</td><td>3.92 × 10<sup>-5</sup></td></tr><tr><td>DNAJC15</td><td>rs17553284</td><td>6.26 × 10<sup>-6</sup></td><td>4.43 × 10<sup>-3</sup></td><td>-19.484</td><td>0.469141</td><td>0.2646227</td><td>0.0112012</td><td>0.0594694</td><td>2.67 × 10<sup>-5</sup></td></tr><tr><td>AGA</td><td>rs4690523</td><td>9.97 × 10<sup>-6</sup></td><td>6.57 × 10<sup>-3</sup></td><td>-16.364</td><td>0.50371342</td><td>0.14900309</td><td>0.06567993</td><td>0.18234469</td><td>7.27 × 10<sup>-5</sup></td></tr><tr><td>TAS2R64P</td><td>rs7297949</td><td>1.05 × 10<sup>-5</sup></td><td>6.85 × 10<sup>-3</sup></td><td>-16.077</td><td>0.0824359</td><td>0.02563944</td><td>0.86790276</td><td>0.878125</td><td>2.53 × 10<sup>-4</sup></td></tr><tr><td>AC009542.2</td><td>rs10954493</td><td>1.07 × 10<sup>-5</sup></td><td>6.95 × 10<sup>-3</sup></td><td>-18.232</td><td>0.97990037</td><td>0.03550566</td><td>0.05021882</td><td>0.21802336</td><td>7.21 × 10<sup>-6</sup></td></tr><tr><td>GTSF1</td><td>rs11170917</td><td>1.46 × 10<sup>-5</sup></td><td>8.84 × 10<sup>-3</sup></td><td>22.336</td><td>0.34080682</td><td>0.02176511</td><td>0.67946143</td><td>0.83764288</td><td>2.02 × 10<sup>-3</sup></td></tr><tr><td>CBLB</td><td>rs1503920</td><td>1.50 × 10<sup>-5</sup></td><td>9.02 × 10<sup>-3</sup></td><td>14.714</td><td>0.02302123</td><td>0.04722184</td><td>0.47428316</td><td>0.98231439</td><td>1.98 × 10<sup>-3</sup></td></tr><tr><td>SOHLH2</td><td>rs9575600</td><td>1.65 × 10<sup>-5</sup></td><td>9.76 × 10<sup>-3</sup></td><td>16.140</td><td>0.10582542</td><td>0.01653465</td><td>0.87063264</td><td>0.92481755</td><td>3.25 × 10<sup>-4</sup></td></tr><tr><td>KANSL1-AS1</td><td>rs17585974</td><td>1.67 × 10<sup>-5</sup></td><td>9.81 × 10<sup>-3</sup></td><td>23.860</td><td>0.11809498</td><td>0.06863156</td><td>0.62204626</td><td>0.77192668</td><td>4.91 × 10<sup>-3</sup></td></tr><tr><td>NSUN7</td><td>rs2437317</td><td>1.73 × 10<sup>-5</sup></td><td>9.94 × 10<sup>-3</sup></td><td>23.457</td><td>0.01246184</td><td>0.01516865</td><td>0.08784461</td><td>0.56754666</td><td>3.00 × 10<sup>-3</sup></td></tr><tr><td>NSF</td><td>rs17698176</td><td>1.74 × 10<sup>-5</sup></td><td>9.97 × 10<sup>-3</sup></td><td>20.367</td><td>0.22639048</td><td>0.00414694</td><td>0.7413922</td><td>0.79027276</td><td>4.09 × 10<sup>-5</sup></td></tr><tr><td>NDUFS5</td><td>rs10888650</td><td>1.82 × 10<sup>-5</sup></td><td>1.01 × 10<sup>-2</sup></td><td>-13.620</td><td>0.52728811</td><td>0.42320815</td><td>0.24332657</td><td>0.03161536</td><td>8.82 × 10<sup>-6</sup></td></tr><tr><td>CASQ1</td><td>rs6693877</td><td>1.98 × 10<sup>-5</sup></td><td>1.10 × 10<sup>-2</sup></td><td>-14.309</td><td>0.1020241</td><td>0.0981154</td><td>0.9077578</td><td>0.3968772</td><td>6.41 × 10<sup>-4</sup></td></tr><tr><td>NPM1P35</td><td>rs4488202</td><td>2.38 × 10<sup>-5</sup></td><td>1.29 × 10<sup>-2</sup></td><td>-21.404</td><td>0.6814926</td><td>0.2694253</td><td>0.1968538</td><td>0.1925592</td><td>1.02 × 10<sup>-3</sup></td></tr><tr><td>PTTG1</td><td>rs1105789</td><td>2.43 × 10<sup>-5</sup></td><td>1.32 × 10<sup>-2</sup></td><td>15.471</td><td>0.0267482</td><td>0.01723734</td><td>0.54238551</td><td>0.79645133</td><td>7.37 × 10<sup>-4</sup></td></tr><tr><td>C8orf46</td><td>rs12541098</td><td>4.07 × 10<sup>-5</sup></td><td>1.98 × 10<sup>-2</sup></td><td>-15.621</td><td>0.01999623</td><td>0.01458403</td><td>0.4383046</td><td>0.67636289</td><td>2.10 × 10<sup>-3</sup></td></tr><tr><td>C8orf89</td><td>rs6472746</td><td>4.65 × 10<sup>-5</sup></td><td>2.19 × 10<sup>-2</sup></td><td>-17.620</td><td>0.004634949</td><td>0.000225028</td><td>0.080493823</td><td>0.022123942</td><td>3.17 × 10<sup>-4</sup></td></tr><tr><td>LINC00654</td><td>rs6116828</td><td>4.93 × 10<sup>-5</sup></td><td>2.29 × 10<sup>-2</sup></td><td>12.740</td><td>0.16132241</td><td>0.02684875</td><td>0.56544041</td><td>0.76958675</td><td>4.19 × 10<sup>-4</sup></td></tr><tr><td>MTPAP</td><td>rs1762598</td><td>5.02 × 10<sup>-5</sup></td><td>2.31 × 10<sup>-2</sup></td><td>-15.032</td><td>0.05145106</td><td>0.01354065</td><td>0.49610048</td><td>0.4799675</td><td>4.82 × 10<sup>-3</sup></td></tr><tr><td>CARD17</td><td>rs1623342</td><td>5.17 × 10<sup>-5</sup></td><td>2.36 × 10<sup>-2</sup></td><td>15.497</td><td>0.30416008</td><td>0.03915456</td><td>0.60188393</td><td>0.84825756</td><td>1.63 × 10<sup>-3</sup></td></tr></tbody></table>
<table><caption>**Table 4.** *Cont.*</caption><thead><tr><th rowspan="2">eGene</th><th rowspan="2">Best-eSNP</th><th colspan="3">From MatrixEQTL Analysis MINUS DISEASE</th><th colspan="5">From Linear Regression to Test Genotype-Disease Interaction</th></tr><tr><th>p-Value</th><th>FDR</th><th>beta</th><th>Intercept p-Val</th><th>eSNP p-Val</th><th>Group p-Val</th><th>eSNP*Group p-Val</th><th>Model p-Val</th></tr></thead><tbody><tr><td>SLC13A5</td><td>rs9889374</td><td>5.49 × 10<sup>-5</sup></td><td>2.48 × 10<sup>-2</sup></td><td>15.823</td><td>0.10733714</td><td>0.53938538</td><td>0.84592422</td><td>0.05661093</td><td>4.53 × 10<sup>-4</sup></td></tr><tr><td>LINC01273</td><td>rs6020255</td><td>5.91 × 10<sup>-5</sup></td><td>2.66 × 10<sup>-2</sup></td><td>-14.134</td><td>0.4831266</td><td>0.461141</td><td>0.2441706</td><td>0.1435129</td><td>1.01 × 10<sup>-3</sup></td></tr><tr><td>LOC101928775</td><td>rs10982832</td><td>5.95 × 10<sup>-5</sup></td><td>2.66 × 10<sup>-2</sup></td><td>16.227</td><td>0.19619467</td><td>0.01188732</td><td>0.8784099</td><td>0.98711196</td><td>1.54 × 10<sup>-4</sup></td></tr><tr><td>VN1R84P</td><td>rs2015481</td><td>6.03 × 10<sup>-5</sup></td><td>2.68 × 10<sup>-2</sup></td><td>13.299</td><td>0.01188751</td><td>0.10670774</td><td>0.3118514</td><td>0.42971681</td><td>2.99 × 10<sup>-4</sup></td></tr><tr><td>CPS1</td><td>rs2250976</td><td>7.35 × 10<sup>-5</sup></td><td>3.09 × 10<sup>-2</sup></td><td>14.858</td><td>0.19187507</td><td>0.09947661</td><td>0.68926108</td><td>0.80062325</td><td>7.05 × 10<sup>-3</sup></td></tr><tr><td>FGF12-AS1</td><td>rs10937543</td><td>8.44 × 10<sup>-5</sup></td><td>3.44 × 10<sup>-2</sup></td><td>-17.045</td><td>0.53243677</td><td>0.07383697</td><td>0.22123134</td><td>0.52782373</td><td>4.57 × 10<sup>-4</sup></td></tr><tr><td>MYRF</td><td>rs7925523</td><td>1.02 × 10<sup>-4</sup></td><td>4.12 × 10<sup>-2</sup></td><td>18.305</td><td>0.049161643</td><td>0.000625781</td><td>0.393915774</td><td>0.080254068</td><td>5.05 × 10<sup>-4</sup></td></tr><tr><td>ITPKB-IT1</td><td>rs697845</td><td>1.03 × 10<sup>-4</sup></td><td>4.15 × 10<sup>-2</sup></td><td>19.480</td><td>0.006181248</td><td>0.000344546</td><td>0.057691247</td><td>0.056088085</td><td>5.11 × 10<sup>-4</sup></td></tr><tr><td>MIR4299</td><td>rs7126296</td><td>1.04 × 10<sup>-4</sup></td><td>4.15 × 10<sup>-2</sup></td><td>15.096</td><td>0.03207948</td><td>0.01109053</td><td>0.62800608</td><td>0.46298706</td><td>3.04 × 10<sup>-3</sup></td></tr><tr><td>CITF22-49E9.3</td><td>rs137878</td><td>1.12 × 10<sup>-4</sup></td><td>4.38 × 10<sup>-2</sup></td><td>-17.668</td><td>0.20391064</td><td>0.02392672</td><td>0.9198819</td><td>0.67454661</td><td>4.21 × 10<sup>-3</sup></td></tr><tr><td>SLC30A8</td><td>rs10505312</td><td>1.18 × 10<sup>-4</sup></td><td>4.57 × 10<sup>-2</sup></td><td>20.222</td><td>0.7540111</td><td>0.457554</td><td>0.3942176</td><td>0.4400821</td><td>3.34 × 10<sup>-2</sup></td></tr><tr><td>CTC-482H14.5</td><td>rs2620833</td><td>1.25 × 10<sup>-4</sup></td><td>4.74 × 10<sup>-2</sup></td><td>16.078</td><td>0.084645173</td><td>0.008758073</td><td>0.852309458</td><td>0.531835913</td><td>1.04 × 10<sup>-3</sup></td></tr><tr><td>PROSER2-AS1</td><td>rs7900122</td><td>1.26 × 10<sup>-4</sup></td><td>4.75 × 10<sup>-2</sup></td><td>-15.424</td><td>0.1660911</td><td>0.03140426</td><td>0.69555964</td><td>0.72078802</td><td>1.49 × 10<sup>-2</sup></td></tr><tr><td>PKN3</td><td>rs10819449</td><td>1.32 × 10<sup>-4</sup></td><td>4.91 × 10<sup>-2</sup></td><td>19.043</td><td>0.3171833</td><td>0.416917</td><td>0.9639612</td><td>0.2775942</td><td>9.58 × 10<sup>-3</sup></td></tr></tbody></table>
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Boxplots of eGene expression in relation to the best-eSNP genotype, stratified by Group: Control vs Disease. eGenes C8orf89, CASQ1, ERICH1. The genotype of each SNP is expressed as 0, 1, 2, and the gene expression data corresponds to the raw data as mean fluorescence value. CTRL = Control; IUGR = Intra-uterine Growth Restriction; PE = Preeclampsia. Genes discussed in the text, because of their relevance to placental biology, are highlighted in green.</figcaption></figure>

**4. Discussion**

In this paper, we attempted a novel approach to uncover relevant associations between genetic variants and gene expression in the placenta from healthy and pathological pregnancies. We started by carefully identifying covariates that may have masked important results in previous studies. As observed previously, gestational age is a major covariate and is strongly associated with disease status. In PCA analysis, the fifth axis is the most strongly associated with gestational age (r = 0.49, Figure 2A), and the most strongly associated to placental disease (r = -0.47, Figure 2A). The inverse correlation between the two features is consistently seen in many transcriptomic studies comparing normal and pathological placentas. The problem is that, almost by definition, placentas from pregnancy complications diseases are usually obtained earlier in gestation than from normal pregnancies. The originality of our approach is as follows: we performed two separate searches for eQTLs, one using all covariates available (including disease status) and one without the disease status used as a covariate, assuming that the difference in the discovery will be due more or less exclusively to the disease. At the initiation of this work, we performed a thorough characterization of the covariates influencing gene expression, and while we could not fully distinguish between “disease group” and “gestational age” given their strong correlation with PC5, we could partially account for gestational age confounding by statistically correcting for PC2 which is strongly correlated (0.39, *p* val < 0.01), but not with “disease group”.

In the first part of our results which included “disease group” as a covariate, we expanded the current knowledge of genetic modulation of gene expression in the placenta, adding to the relatively limited number of previous studies in the field [7–9,12]. We identified 279 placental *cis*-eQTLs (FDR ≤ 0.05) in 57 human placental samples, correlating with expression changes of 43 unique eGenes, among which 22 constituted novel eGenes over the previous studies. One hundred and eleven (111) of the total 279 *cis*-eQTLs, including the most significant eSNPs, were associated with previously described placental eGenes ZSCAN9, PSG7, TOB2P1, ERAP2, and AQP11, corroborating the validity of our analysis. Interestingly, for the first time in the placenta, we describe *cis*-eQTL eSNPs that correlate with neighboring transcripts, identifying potential placental regulatory hotspots. On chromosome 11, a novel eGene LOC646029 and downstream AQP11 shared 10 eSNPs, while the three eSNPs associated with CARD17 expression levels were also eSNPs for novel placental eGene CASP1P2. Similarly, 24 out of the 42 eSNPs associated with ZSCAN9
were found to correlate with changes in expression of the upstream transcript *TOB2P1* on chromosome 6. *CASP1P2* and *TOB2P1* are both pseudogenes, which implies that they may not encompass functional open reading frames and that changes in their expression could not have an impact on placental function as proteins. However, the changes in their expression levels as well as in neighboring genes as a result of genetic changes, highlight these regions as potential placental regulatory sites. As additional evidence of the relevance of these regions to placental function, genotype changes in these loci show a statistically significant interaction with the disease status, resulting in significant changes in the expression of ZSCAN9 and LOC646029, as well as AQP11.

Kikas and coworkers extensively reviewed previous studies on genetic regulation of gene expression in the placenta, including gene-candidate approaches, as well as the three genome-wide studies mentioned above [40]. Peng and coworkers combined RNA-sequencing and genotype data of 159 healthy placentas, laying the foundations for eQTL analyses in the placenta [9]. They report only the top eSNP of each eGene, therefore identifying 3218 *cis*-eQTLs correlating with 3218 unique transcripts and 35 *trans*-eQTLs (*cis*-distance between eSNP and eGene < 500kb) with an FDR threshold, defined by permutations, of 0.1%. This limited number of *trans*-QTLs, even with a sample size that is more than double one of other studies in placental eQTLs, prompted to restrict further efforts to *cis*-QTL identification alone. Delahaye and coworkers, expanded the repertoire of placental *cis*-eQTL by identifying 985 eSNPs, corresponding to 615 eGenes in 80 healthy placentas, as well as investigating methylation-QTLs in 300 placenta samples [7]. Finally, Kikas and coworkers, with a total of 40 samples, identified 199 significant *cis*-eQTLs (FDR < 0.05) corresponding to 63 unique genes in a window of 100 kb around each gene [8]. Both studies showed considerable overlap with [9], respectively, for 62% and 80% of the eGenes. However, if we consider the percentage of unique eGenes, not replicated across all studies, 88%, 37%, and 21% of eGenes were identified, respectively, by Peng et al. [9], Delahaye et al. [7] and Kikas et al. [8].

When comparing our results with previous studies, the overall level of overlap between eGene lists is 51%, which could be due to differences in population stratification, sample size as well as differences in analysis workflow. In particular, the number of overlapping eGenes is n = seventeen for the [9] dataset, n = eight for [7], and n = five for [8]. The most significant overlap is with the work of [8]; Fisher's exact tests $p < 2.75 \times 10^{-8}$, 1-2 orders of magnitude more significant than the overlap with the two other studies. This could be explained by the fact that in terms of analytical design, our study is closer to the one of [8]. In both analyses, cis-QTLs were analyzed in a window of 100 kb, for a total number of variants of 417,114, in the same order of magnitude, similarly; the sample sizes are also close with 40 and 57 placentas in [8] and our study respectively. Probably the most crucial similarities are the use of the same software for the eQTL analysis (i.e., MatrixEQTL in R [24]) and the inclusion of placentas from pregnancies complicated by PE. The use of the same analytical software improves the replicability of the results since the statistical tests are the same and can have a huge impact on which eSNP-eGene pairs pass the threshold of significance.

Similarly, the presence of PE samples implies a similar pattern of gene expression; however, in our study design, we did consider "disease group" as a confounding factor and included it as a covariate, which should reduce the effects on the output.

With our analysis comparing the two datasets, the ALL COVARIATES vs. MINUS DISEASE, that differ only by removing the disease as a confounding covariate in the second, we propose a new strategy for finding eQTL that are relevant in pathological pregnancies. In particular, this approach allowed us first to validate our study design as well as sample quality by being able to compare our findings with the existing literature. In the second instance, we were able to further investigate the role of genotype changes on the expression of genes for which disease has a relevant influence and could therefore be important for placental function.
Genome-Wide Association Studies (GWAS) use large sample sizes (in the tens of hundreds) in an attempt to identify the enrichment of variant alleles in the disease group compared to healthy controls. GWAS conducted in the placenta has identified SNPs upstream of the FLT1 locus, rs4769613, as harboring the risk allele (T) associated with PE [41]. Despite the limited functional understanding offered by GWAS approaches, the use of combinatorial approaches that exploit knowledge from different types of datasets helps draw these functional conclusions [42–44]. Recently, the rs4769613 locus was genotyped in a cohort of placental samples, including 57 PE and 277 control; the presence of the risk allele correlated with changes in the placental gene expression level of FLT1 only in the PE group, defining rs4769613 as a conditional eQTL [31]. This suggests that the risk variant, when present in the placentas of mothers with PE, favors the increase in expression of the FLT1 gene and, therefore, could be located in a regulatory region for the gene itself, targeted by regulatory factors that are expressed specifically in the PE placenta.

In our work, we identified multiple genes that presented a promising 'genotype and disease' interaction, which could therefore represent placental conditional eQTLs relevant to disease etiology. The genes found during this stage are supposed to be influenced by the disease either exclusively or mainly but with an influence of other factors.

Among these, AQP11, IL36R, FUT10, PTTG1, CBLB, and DNAJC15 are interesting both in terms of behavior, as well as the documented role in placental function.

Among the 'disease exclusive' genes represented in detail (in green boxes), there are PTTG1, DNAJC15 (Figure 6), TAS2R64P and LINC00654 (Figure 7), C8ORF89 and ERICH1 (Figure 8). The profiles of expression are very different; for instance, in the case of ERICH1, the SNP influence of the profile is the same way in disease and controls, but overall in the disease, the expression is systematically inferior. By contrast, C8ORF89 has an SNP-dependent profile in controls but not in disease; PTTG1 is overexpressed in disease but only rather in specific SNP configuration.

PTTG1 has been involved in trophoblast invasion through the regulation of matrix metalloproteases MMP2 and MMP9 by regulating integrin/rho signalling [34]. AQP11 belongs to the family of water transporters through the cellular and endoplasmic reticulum membrane. It has recently been implicated in redox regulation, and other members of the same family, such as AQP, were found to be involved in the mechanism of aspirin action in preventing PE [45–47]. IL36RN is highly expressed in the placenta, where it could be involved in trophoblast proliferation [48]. FUT10, involved in the fucosylation of proteins, has been shown to promote the binding of TNF to its receptor, thus promoting the inflammatory reaction, which is relevant for placental diseases [49]. Interestingly, we identified a significant interaction between disease and genotype on the levels of FUT10. DNAJC15 is reported as an important factor in mitochondrial function and regulation of the respiratory chain. It has been abundantly demonstrated by us and many others how these cascades of oxidative phosphorylation and oxidative stress are pivotal in placental diseases [50,51]. CBLB has been recently shown to influence Natural Killer cell differentiation [52] and Uterine NK cells are pivotal for efficient implantation. In summary, all these genes may be involved in placental disease by a contribution to various pathways, which explains why they are not explicitly regrouped in a biological function.

**5. Study Limitations and Conclusions**

Amongst the limitations of our study, there is, of course, the sample size, which makes our data a preliminary incitement that needs confirmation in further studies. The choice to subtract extensively the covariate effect encompasses a risk of missing positive hits (Type II error), but our objective was to decrease the Type I error and minimize the finding of false positives. Some additional analyses could be performed by removing some of the covariates to detect eSNPs influenced by sex, but this is beyond the objective of our current work. Furthermore, among the control placentas, one-third were from natural deliveries, which we included in our study by adding the delivery mode as an additional covariate). Gestational age was found to induce 10% of the variability, consistent with
previous observations in the field. Subtracting its effect and using it as a covariate might be considered a mathematical trick that could mask complex interactions between the disease effects and the gestational age effects. Nevertheless, it is currently difficult to adopt a better strategy in the field of placental disease.

The approach based on two-hits, using unusually the disease as a covariate and then adding it afterward, allows by subtraction to identify novel couples of SNP-eGenes that may be of interest for placental disease studies or other cases of human diseases where confounding factors are important. We hope that our work will be a starting point for researchers using other complex datasets available for analysis.

**Supplementary Materials:** The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/cells12040578/s1, Figure S1, Results of the Principal Component Analysis performed to assess population stratification and establish admixture. Table S1 present the parameters for threshold evaluation of the residual variance, Table S2 give the dataset that was used for the analyses, Tables S3 and S4 list statistically significant couples eSNP-eGene for the dataset including disease as a covariate and for the one excluding it. Supplementary Tables S5 and S6 lists the input genes with their expression levels, and residual variances, respectively. Table S7 ists the top-eSNP in terms of statistical significance.

**Author Contributions:** Conceptualization, D.V., C.E.M., B.T. and A.K.; methodology, D.V., C.A., F.M., C.M., S.J., C.G.; software, C.A., C.S.M.R., F.M.; validation, C.A., D.V.; formal analysis, D.V., C.A.; investigation, C.A.; resources, B.T., A.K., V.G., G.G., L.M., C.M.; data curation, C.A., D.V.; writing—original draft preparation, C.A., D.V.; writing—review and editing, C.A., B.T., A.K., C.E.M., F.M., C.M., D.V.; visualization, D.V.; supervision, D.V.; project administration, C.M., D.V.; funding acquisition, C.E.M., B.T., A.K., G.G., C.M., D.V. All authors have read and agreed to the published version of the manuscript.

**Funding:** This project was funded by European Union's Horizon 2020 research and innovation program under Marie Skłodowska-Curie Actions Innovative Training Network (H2020-MSCA-ITN 2017), Grant No. 765274, acronym IPLACENTA (http://www.iplacenta.eu, accessed on 27 January 2023).

**Institutional Review Board Statement:** For the London StGeorge's cohort, the administrative ethical references were 19/LO/0974, approved by the Brent Research Ethic Committee (London) on 28 June 2019. For the Angers cohort, the study was approved by the Ethics Committee of Angers. The cohort was registered at the French CNIL (Commission Nationale de l'Informatique et des Libertés no. pWP03752UL, ethics committee for the collection of clinical data from patient records). The study was validated by the French CPP (Comité de Protection des Personnes) and registered to the French Ministry of Research under number DC-2009-907. For the Cochin hospital cohort, the ethical administration was given by the CPP 'Ile de France XI' under the reference number 11018, 3 March 2011.

**Informed Consent Statement:** Informed consent was obtained from all subjects involved in the study.

**Data Availability Statement:** The significant data of the study are provided in detail in the Supplementary Tables. Access to the complete dataset can be obtained on request from the authors.

**Acknowledgments:** We thank the genomic platform of Cochin Institute (Paris) and the Post-Genomic Platform of Sorbonne University (Paris), where the genomic data were generated. We thank Odile Blanchet and the Centre for Biological Resources (CHU Angers, France) for the provision of biological samples from the Angers University Hospital Cohort. We would also like to thank Etienne Patin, Aurélie Cobat and Maxime Rotival for their valuable insights and helpful advice.

**Conflicts of Interest:** The authors declare no conflict of interest.

## References

1. Griffith, O.W.; Chavan, A.R.; Protopapas, S.; Maziarz, J.; Romero, R.; Wagner, G.P. Embryo implantation evolved from an ancestral inflammatory attachment reaction. *Proc. Natl. Acad. Sci. USA* **2017**, *114*, E6566–E6575. [CrossRef] [PubMed]

2. Burton, G.J.; Fowden, A.L.; Thornburg, K.L. Placental Origins of Chronic Disease. *Physiol. Rev.* **2016**, *96*, 1509–1565. [CrossRef] [PubMed]
3. Desoye, G.; Hauguel-de Mouzon, S. The human placenta in gestational diabetes mellitus. The insulin and cytokine network. *Diabetes Care* **2007**, *30* (Suppl. S2), S120-S126. [CrossRef] [PubMed]

4. Lonsdale, J.; Thomas, J.; Salvatore, M.; Phillips, R.; Lo, E.; Shad, S.; Hasz, R.; Walters, G.; Garcia, F.; Young, N.; et al. The Genotype-Tissue Expression (GTEx) project. *Nat. Genet.* **2013**, *45*, 580–585. [CrossRef]

5. Leavey, K.; Benton, S.J.; Grynspan, D.; Kingdom, J.C.; Bainbridge, S.A.; Cox, B.J. Unsupervised Placental Gene Expression Profiling Identifies Clinically Relevant Subclasses of Human Preeclampsia. *Hypertension* **2016**, *68*, 137–147. [CrossRef]

6. Madeleneau, D.; Buffat, C.; Mondon, F.; Grimault, H.; Rigourd, V.; Tsatsaris, V.; Letourneur, F.; Vaiman, D.; Barbaux, S.; Gascoin, G. Transcriptomic analysis of human placenta in intrauterine growth restriction. *Pediatr. Res.* **2015**, *77*, 799–807. [CrossRef]

7. Delahaye, F.; Do, C.; Kong, Y.; Ashkar, R.; Salas, M.; Tycko, B.; Wapner, R.; Hughes, F. Genetic variants influence on the placenta regulatory landscape. *PLoS Genet.* **2018**, *14*, e1007785. [CrossRef]

8. Kikas, T.; Rull, K.; Beaumont, R.N.; Freathy, R.M.; Laan, M. The Effect of Genetic Variation on the Placental Transcriptome in Humans. *Front. Genet.* **2019**, *10*, 550. [CrossRef]

9. Peng, S.; Deyssenroth, M.A.; Di Narzo, A.F.; Lambertini, L.; Marsit, C.J.; Chen, J.; Hao, K. Expression quantitative trait loci (eQTLs) in human placentas suggest developmental origins of complex diseases. *Hum. Mol. Genet.* **2017**, *26*, 3432–3441. [CrossRef]

10. Leek, J.T.; Johnson, W.E.; Parker, H.S.; Jaffe, A.E.; Storey, J.D. The sva package for removing batch effects and other unwanted variation in high-throughput experiments. *Bioinformatics* **2012**, *28*, 882–883. [CrossRef]

11. Bergstedt, J.; Azzou, S.A.K.; Tsuo, K.; Jaquaniello, A.; Urrutia, A.; Rotival, M.; Lin, D.T.S.; MacIsaac, J.L.; Kobor, M.S.; Albert, M.L.; et al. The immune factors driving DNA methylation variation in human blood. *Nat. Commun.* **2022**, *13*, 5895. [CrossRef]

12. Ruano, C.S.M.; Apicella, C.; Jacques, S.; Gascoin, G.; Gaspar, C.; Miralles, F.; Méhats, C.; Vaiman, D. Alternative splicing in normal and pathological human placentas is correlated to genetic variants. *Hum. Genet.* **2021**, *140*, 827–848. [CrossRef]

13. Chabrun, F.; Huetz, N.; Dieu, X.; Rousseau, G.; Bouzille, G.; Chao de la Barca, J.M.; Procaccio, V.; Lenaers, G.; Blanchet, O.; Legendre, G.; et al. Data-Mining Approach on Transcriptomics and Methylomics Placental Analysis Highlights Genes in Fetal Growth Restriction. *Front. Genet.* **2019**, *10*, 1292. [CrossRef]

14. Gascoin-Lachambre, G.; Buffat, C.; Rebourcet, R.; Chelbi, S.T.; Rigourd, V.; Mondon, F.; Mignot, T.M.; Legras, E.; Simeoni, U.; Vaiman, D.; et al. Cullins in human intra-uterine growth restriction: Expressional and epigenetic alterations. *Placenta* **2010**, *31*, 151–157. [CrossRef]

15. Brown, M.A.; Magee, L.A.; Kenny, L.C.; Karumanchi, S.A.; McCarthy, F.P.; Saito, S.; Hall, D.R.; Warren, C.E.; Adoyi, G.; Ishaku, S.; et al. The hypertensive disorders of pregnancy: ISSHP classification, diagnosis & management recommendations for international practice. *Pregnancy Hypertens.* **2018**, *13*, 291–310. [CrossRef]

16. Irizarry, R.A.; Hobbs, B.; Collin, F.; Beazer-Barclay, Y.D.; Antonellis, K.J.; Scherf, U.; Speed, T.P. Exploration, normalization, and summaries of high density oligonucleotide array probe level data. *Biostatistics* **2003**, *4*, 249–264. [CrossRef]

17. Lupattelli, A.; Wood, M.E.; Nordeng, H. Analyzing Missing Data in Perinatal Pharmacoepidemiology Research: Methodological Considerations to Limit the Risk of Bias. *Clin. Ther.* **2019**, *41*, 2477–2487. [CrossRef]

18. Blighe, K.; Lun, A. PCAtools: PCAtools: Everything Principal Components Analysis [WWW Document]. R Packag. Version 2.4.0. 2021. Available online: https://github.com/kevinblighe/PCAtools (accessed on 9 January 2021).

19. Tabachnick, B.G.; Fidell, L.S. *Using Multivariate Statistics*, 7th ed.; Pearson: Boston, MA, USA, 2019.

20. Chang, C.C.; Chow, C.C.; Tellier, L.C.; Vattikuti, S.; Purcell, S.M.; Lee, J.J. Second-generation PLINK: Rising to the challenge of larger and richer datasets. *Gigascience* **2015**, *4*, 7. [CrossRef]

21. Purcell, S.; Neale, B.; Todd-Brown, K.; Thomas, L.; Ferreira, M.A.; Bender, D.; Maller, J.; Sklar, P.; De Bakker, P.I.; Daly, M.J.; et al. PLINK: A tool set for whole-genome association and population-based linkage analyses. *Am. J. Hum. Genet.* **2007**, *81*, 559–575. [CrossRef]

22. Auton, A.; Brooks, L.D.; Durbin, R.M.; Garrison, E.P.; Kang, H.M.; Korbel, J.O.; Marchini, J.L.; McCarthy, S.; McVean, G.A.; Abecasis, G.R. A global reference for human genetic variation. *Nature* **2015**, *526*, 68–74. [CrossRef]

23. Sudmant, P.H.; Rausch, T.; Gardner, E.J.; Handsaker, R.E.; Abyzov, A.; Huddleston, J.; Zhang, Y.; Ye, K.; Jun, G.; Fritz, M.H.; et al. An integrated map of structural variation in 2,504 human genomes. *Nature* **2015**, *526*, 75–81. [CrossRef] [PubMed]

24. Shabalin, A.A. Matrix eQTL: Ultra fast eQTL analysis via large matrix operations. *Bioinformatics* **2012**, *28*, 1353–1358. [CrossRef] [PubMed]

25. Listgarten, J.; Kadie, C.; Schadt, E.E.; Heckerman, D. Correction for hidden confounders in the genetic analysis of gene expression. *Proc. Natl. Acad. Sci. USA* **2010**, *107*, 16465–16470. [CrossRef] [PubMed]

26. Mele, M.; Ferreira, P.G.; Reverter, F.; DeLuca, D.S.; Monlong, J.; Sammeth, M.; Young, T.R.; Goldmann, J.M.; Pervouchine, D.D.; Sullivan, T.J.; et al. Human genomics. The human transcriptome across tissues and individuals. *Science* **2015**, *348*, 660–665. [CrossRef] [PubMed]

27. Rouhani, F.; Kumasaka, N.; De Brito, M.C.; Bradley, A.; Vallier, L.; Gaffney, D. Genetic background drives transcriptional variation in human induced pluripotent stem cells. *PLoS Genet.* **2014**, *10*, e1004432. [CrossRef]

28. Rivals, I.; Personnaz, L.; Taing, L.; Potier, M.C. Enrichment or depletion of a GO category within a class of genes: Which test? *Bioinformatics* **2007**, *23*, 401–407. [CrossRef]

29. Gong, S.; Gaccioli, F.; Dopierała, J.; Sovio, U.; Cook, E.; Volders, P.J.; Martens, L.; Kirk, P.D.W.; Richardson, S.; Smith, G.C.S.; et al. The RNA landscape of the human placenta in health and disease. *Nat. Commun.* **2021**, *12*, 2639. [CrossRef]

30. Willyard, C. New human gene tally reignites debate. *Nature* **2018**, *558*, 354–355. [CrossRef]
31. Kikas, T.; Inno, R.; Ratnik, K.; Rull, K.; Laan, M. C-allele of rs4769613 Near FLT1 Represents a High-Confidence Placental Risk Factor for Preeclampsia. *Hypertension* **2020**, *76*, 884–891. [CrossRef]

32. Yan, Y.H.; Yi, P.; Zheng, Y.R.; Yu, L.L.; Han, J.; Han, X.M.; Li, L. Screening for preeclampsia pathogenesis related genes. *Eur. Rev. Med. Pharmacol. Sci.* **2013**, *17*, 3083–3094.

33. May, A.P.; Whiteheart, S.W.; Weis, W.I. Unraveling the mechanism of the vesicle transport ATPase NSF, the N-ethylmaleimide-sensitive factor. *J. Biol. Chem.* **2001**, *276*, 21991–21994. [CrossRef]

34. Lim, S.M.; Jang, H.Y.; Lee, J.E.; Shin, J.S.; Park, S.H.; Yoon, B.H.; Kim, G.J. Alteration of Pituitary Tumor Transforming Gene-1 Regulates Trophoblast Invasion via the Integrin/Rho-Family Signaling Pathway. *PLoS ONE* **2016**, *11*, e0149371. [CrossRef]

35. Liu, R.; Wang, X.; Shen, Y.; He, A. Long non-coding RNA-based glycolysis-targeted cancer therapy: Feasibility, progression and limitations. *Mol. Biol. Rep.* **2021**, *48*, 2713–2727. [CrossRef]

36. Steinman, G. The putative etiology and prevention of autism. *Prog. Mol. Biol. Transl. Sci.* **2020**, *173*, 1–34. [CrossRef]

37. Stuebe, A.M.; Lyon, H.; Herring, A.H.; Ghosh, J.; Wise, A.; North, K.E.; Siega-Riz, A.M. Obesity and diabetes genetic variants associated with gestational weight gain. *Am. J. Obstet. Gynecol.* **2010**, *203*, 283.e1–283.e17. [CrossRef]

38. Hamada, H.; Okae, H.; Toh, H.; Chiba, H.; Hiura, H.; Shirane, K.; Sato, T.; Suyama, M.; Yaegashi, N.; Sasaki, H.; et al. Allele-Specific Methylome and Transcriptome Analysis Reveals Widespread Imprinting in the Human Placenta. *Am. J. Hum. Genet.* **2016**, *99*, 1045–1058. [CrossRef]

39. Möpert, K.; Löffler, K.; Röder, N.; Kaufmann, J.; Santel, A. Depletion of protein kinase N3 (PKN3) impairs actin and adherens junctions dynamics and attenuates endothelial cell activation. *Eur. J. Cell Biol.* **2012**, *91*, 694–705. [CrossRef]

40. Kikas, T.; Laan, M.; Kasak, L. Current knowledge on genetic variants shaping placental transcriptome and their link to gestational and postnatal health. *Placenta* **2021**, *116*, 2–11. [CrossRef]

41. McGinnis, R.; Steinthorsdottir, V.; Williams, N.O.; Thorleifsson, G.; Shooter, S.; Hjartardottir, S.; Bumpstead, S.; Stefansdottir, L.; Hildyard, L.; Sigurdsson, J.K.; et al. Variants in the fetal genome near FLT1 are associated with risk of preeclampsia. *Nat. Genet.* **2017**, *49*, 1255–1260. [CrossRef]

42. Ishigaki, K. Beyond GWAS: From simple associations to functional insights. *Semin. Immunopathol.* **2022**, *44*, 3–14. [CrossRef]

43. Laivuori, H. Pitfalls in setting up genetic studies on preeclampsia. *Pregnancy Hypertens.* **2013**, *3*, 60. [CrossRef] [PubMed]

44. Tam, V.; Patel, N.; Turcotte, M.; Bosse, Y.; Pare, G.; Meyre, D. Benefits and limitations of genome-wide association studies. *Nat. Rev. Genet.* **2019**, *20*, 467–484. [CrossRef] [PubMed]

45. Bestetti, S.; Galli, M.; Sorrentino, I.; Pinton, P.; Rimessi, A.; Sitia, R.; Medrano-Fernandez, I. Human aquaporin-11 guarantees efficient transport of H₂O₂ across the endoplasmic reticulum membrane. *Redox Biol.* **2020**, *28*, 101326. [CrossRef] [PubMed]

46. Hwang, H.S. A new water channel aquaporin-11: Extension to renal transplantation. *Kidney Res. Clin. Pract.* **2015**, *34*, 188–189. [CrossRef]

47. Zhang, J.; Jing, S.; Zhang, H.; Zhang, J.; Xie, H.; Feng, L. Low-dose aspirin prevents LPS-induced preeclampsia-like phenotype via AQP-1 and the MAPK/ERK 1/2 pathway. *Placenta* **2022**, *121*, 61–69. [CrossRef]

48. Lv, Z.; Fan, J.; Zhang, X.; Huang, Q.; Han, J.; Wu, F.; Hu, G.; Guo, M.; Jin, Y. Integrative genomic analysis of interleukin-36RN and its prognostic value in cancer. *Mol. Med. Rep.* **2016**, *13*, 1404–1412. [CrossRef]

49. Yu, H.; Li, M.; Wen, X.; Yang, J.; Liang, X.; Li, X.; Bao, X.; Shu, J.; Ren, X.; Chen, W.; et al. Elevation of alpha-1,3 fucosylation promotes the binding ability of TNFR1 to TNF-alpha and contributes to osteoarthritic cartilage destruction and apoptosis. *Arthritis Res. Ther.* **2022**, *24*, 93. [CrossRef]

50. Aouache, R.; Biquard, L.; Vaiman, D.; Miralles, F. Oxidative Stress in Preeclampsia and Placental Diseases. *Int. J. Mol. Sci.* **2018**, *19*, 1496. [CrossRef]

51. Chatre, L.; Ducat, A.; Spradley, F.T.; Palei, A.C.; Chereau, C.; Couderc, B.; Thomas, K.C.; Wilson, A.R.; Amaral, L.M.; Gaillard, I.; et al. Increased NOS coupling by the metabolite tetrahydrobiopterin (BH4) reduces preeclampsia/IUGR consequences. *Redox Biol.* **2022**, *55*, 102406. [CrossRef]

52. Guo, X.; Mahlakoiv, T.; Ye, Q.; Somanchi, S.; He, S.; Rana, H.; DiFiglia, A.; Gleason, J.; Van der Touw, W.; Hariri, R.; et al. CBLB ablation with CRISPR/Cas9 enhances cytotoxicity of human placental stem cell-derived NK cells for cancer immunotherapy. *J. Immunother. Cancer* **2021**, *9*, e001975. [CrossRef]

**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.