Integrative analysis provides multi‐omics evidence for the pathogenesis of placenta percreta

Abstract Pernicious placenta previa with placenta percreta (PP) is a catastrophic condition during pregnancy. However, the underlying pathogenesis remains unclear. In the present study, the placental tissues of normal cases and PP tissues of pernicious placenta previa cases were collected to determine the expression profile of protein‐coding genes, miRNAs, and lncRNAs through sequencing. Weighted gene co‐expression network analysis (WGCNA), accompanied by miRNA target prediction and correlation analysis, were employed to select potential hub protein‐coding genes and lncRNAs. The expression levels of selected protein‐coding genes, Wnt5A and MAPK13, were determined by quantitative PCR and immunohistochemical staining, and lncRNA PTCHD1‐AS and PAPPA‐AS1 expression levels were determined by quantitative PCR and fluorescence in situ hybridization. The results indicated that 790 protein‐coding genes, 382 miRNAs, and 541 lncRNAs were dysregulated in PP tissues, compared with normal tissues. WGCNA identified coding genes in the module (ME) black and ME turquoise modules that may be involved in the pathogenesis of PP. The selected potential hub protein‐coding genes, Wnt5A and MAPK13, were down‐regulated in PP tissues, and their expression levels were positively correlated with the expression levels of PTCHD1‐AS and PAPPA‐AS1. Further analysis demonstrated that PTCHD1‐AS and PAPPA‐AS1 regulated Wnt5A and MAPK13 expression by interacting with specific miRNAs. Collectively, our results provided multi‐omics data to better understand the pathogenesis of PP and help identify predictive biomarkers and therapeutic targets for PP.


| INTRODUC TI ON
Placenta accreta is an abnormal placental attachment caused by the invasion of placental villi into the myometrium. According to the depth of placental invasion into the myometrium and degree of infiltration into the organs adjacent to the uterus, abnormal placental attachment can occur as follows: (a) placenta accreta: placenta invades the superficial myometrium of the uterus; (b) placenta increta: placenta invades the deep myometrium of the uterus; and (c) placenta percreta (PP): placenta penetrates the uterine wall and reaches the serous layer of the uterus and even invades the organs adjacent to the uterus. 1 The main risk factors for placenta accreta include placenta previa, history of previous caesarean section, history of intrauterine surgery, pregnancy by in vitro fertilization-embryo transfer, advanced age, and history of uterine perforation. [2][3][4] In 1993, Chattopadhyay et al investigated the relationship between placenta previa and placenta accreta and previous caesarean section and put forward the concept of pernicious placenta previa for the first time. 5 Pernicious placenta previa with placenta accreta is an important cause of perinatal hysterectomy, premature delivery, and perinatal death. The hysterectomy rate of patients is as high as 66%, often accompanied by bladder and ureteral injury, and maternal mortality caused by massive haemorrhage is as high as 7%. [6][7][8] The process of placenta accreta development is complex, and its pathogenesis is gaining more attention. Loss of decidua, 9,10 enhanced invasiveness of trophoblasts [11][12][13][14] and abnormal recasting of uterine spiral arteries are considered to be the three important pathophysiological bases that lead to placenta accrete by interacting with and influencing each other. [14][15][16][17] Recently, it has been shown that some important coding and non-coding genes are closely related to placenta accreta, such as those that encode tumour necrosis factor-related apoptosis-inducing ligand-receptor 2 (TRAIL-R2), miR-29 and miR-519d. [18][19][20] However, the regulatory mechanisms of the molecular networks related to placenta accreta remain unclear.
Long non-coding RNAs (lncRNAs) are a class of transcripts longer than 200 base pairs that generally do not code for proteins. LncRNAs exert diverse roles in cellular and biological processes via the regulation of gene expression and chromatin dynamics. 21 Recently, lncRNAs have been shown to contribute to the pathogenesis of various diseases, including cancer, [22][23][24] cardiovascular disease 25 and nervous system diseases. 26 The deregulation of lncRNAs IGF2/H19, MEG3, SPRY4-IT1, HOTAIR, MALAT1, FLT1P1 and CEACAMP8 in placental trophoblasts is involved in the pathogenesis of preeclampsia. 27 Among them, lncRNA MALAT1 may be involved in the pathogenesis of preeclampsia via the regulation of the proliferation, cell cycle, apoptosis, migration and invasion of trophoblast cells. 28 Thus, we speculated that lncRNAs may also be deregulated in the placental tissues of pregnant women with placenta accreta spectrum (PAS) and play a crucial role in the pathogenesis of PAS. Therefore, in this study, the different pathological features of normal placental tissue and penetrating placental tissue were compared. Additionally, the expression profiles of coding genes, lncRNAs and miRNAs in these two types of placental tissue were compared and analysed using multi-omics. Weight gene co-expression network analysis (WGCNA), accompanied by miRNA target prediction and correlation analysis, was employed to select potential hub coding genes and lncRNAs. The expression of the selected coding genes, Wnt5A and MAPK13, was determined by quantitative PCR (qPCR) and immunohistochemical staining and that of lncRNAs PTCHD1-AS and PAPPA-AS1 was determined by qPCR and fluorescence in situ hybridization (FISH) staining. The results have demonstrated the expression profiles of ln-cRNAs, miRNAs and coding genes in PP and provided multi-omics evidence to explain the pathogenesis of placenta accreta.

| RNA extraction and sequencing
Total RNA samples were extracted following the instruction of Trizol kit (Invitrogen) and their quality and quantity were determined as K E Y W O R D S lncRNA, miRNA, pernicious placenta previa, placenta percreta, Wnt5A follows. The samples were first qualified using 1% agarose gel electrophoresis to detect possible contamination and degradation. RNA purity and concentration were then determined using a NanoPhotometer ® spectrophotometer (Implen, Munich, Germany).
Finally, RNA integrity and quantity were measured using the RNA Nano 6000 Assay Kit and the Bioanalyzer 2100 system (Agilent, Santa Clara, CA, USA). RNA sequencing was performed by Chengdu Basebiotech Co., Ltd (Chengdu, China) as follow: a total of 1 μg of RNA per sample was used as the input material for RNA sample preparations. Sequencing libraries were generated using NEBNext ® UltraTM RNA Library Prep Kit for Illumina ® (NEB, Ipswich, MA, USA) following the manufacturer's recommendations, and index codes were added to attribute sequences to each sample. Briefly, mRNA was purified from total RNA using poly T oligo-coupled mag- Finally, PCR products were purified using the AMPure XP system and library quality was assessed on an Agilent Bioanalyzer 2100 system. The clustering of the index-coded samples was performed on a cBot-Cluster Generation System using TruSeq PE Cluster Kit v3-cBot-HS (Illumia) according to the manufacturer's instructions.
After cluster generation, the library preparations were sequenced on an Illumina Novaseq platform and 150 bp paired-end reads were generated. Significantly differentially expressed lncRNAs, miRNAs, and protein-coding genes were screened based on absolute value of log2 (fold change) ≥ 1, at a P value <0.05.

| Weighted gene co-expression network analysis
Using standard weighted gene co-expression network analysis (WGCNA) procedures, a network was constructed using the WGCNA package in R (https://horva th.genet ics.ucla.edu/html/ Coexp ressi onNet work/Rpack ages/WGCNA/) and data were visualized using the Cytoscape software. The flashCluster package in R software was first used to analyse the samples and identify abnormal values. The WGCNA adjacency function was then used to create an adjacency matrix and calculate Pearson's correlations, to determine the consistency of gene expression levels between each gene pair. 29 Next, we used the topological overlap matrix (TOM) similarity function to transform the matrix into a TOM. Finally, co-expression modules were constructed using the WGCNA algorithm, and the gene information of each module was extracted.

| Kyoto Encyclopedia of Genes and Genomes pathway analysis
The Kyoto Encyclopedia of Genes and Genomes (KEGG) is a set of databases that provides a comprehensive understanding of biological systems. It can be used to analyse biological pathways and genes related to diseases and drugs. KEGG pathway datasets were
For lncRNA and protein-coding genes, the default target relationship in mirwalk2 was used.

| Predicting lncRNAs of sponge regulatory network and co-expression module genes
We obtained a list of genes from the clinical trait-related co-expression module. We focused on genes that showed differential expression between normal and diseased tissue. We devised a computational strategy to identify candidate lncRNA-gene pairs based on sponge regulatory network. Firstly, for each lncRNA-driver gene pair, we estimated the significance of shared miRNAs with the same seeds (P-value of one-tailed Fisher's exact test) and the significance of expression correlation across all samples. We then computed a combined P-value by converting the P-values of these two tests, P1 and P2, using the sum of logs method (also called Fisher's method) with the metap package. The candidate lncRNA-RNA-driver gene pairs met the criterion that the adjusted combined P-value was no larger than a threshold of 5% (ie false discovery rate, r < 0.05).
Secondly, we selected lncRNA-protein-coding gene pairs that shared at least ten different miRNAs. Finally, we selected RNA-driver gene pairs that showed at least a moderate positive correlation of their expression levels (r > 0.25). All lncRNAs, genes and miRNAs were filtered to show differential expression at a cutoff P-value < 0.05.
The top 20 miRNAs ranked by P-value were plotted. All of the analysis were conducted according to method report by previous study, 30 with in-house R scripts. ceRNA networks were visualized with Cytoscapesoftwares.

| Venn analysis
To complete the Venn map, we prepared a list of differential genes in each group. We then used the mapping website (http://bioin fogp. cnb.csic.es/tools/ venny/ index.html) to obtain a Venn graph.

| Statistical analysis
All data were analysed by a paired t test, using GraphPad Prism 5.0 (GraphPad software, Inc, San Diego, CA, USA). A P-value < 0.05 was considered statistically significant. All experiments were performed on three or more independent occasions, and the data are presented as the mean ± standard error.

| Clinical characteristics of the collected samples
To investigate the pathogenesis of placenta accreta, placental tissues from pregnant women with normal placenta and those with pernicious placenta previa with PP who met the surgical indications were collected for further sequencing. The detailed clinical characteristics of the pregnant women included in the study are listed in Table 1. Magnetic resonance imaging (MRI) showed that the placenta was located in the normal position of the posterior wall of the uterus without adhesion implantation in the normal group ( Figure 1A). In the PP group, the placenta completely covered the inner cervix, accompanied by abundant blood vessels in the cervix and the disappearance of normal inner cervix morphology ( Figure 1A). The placenta above the inner mouth was 8 cm thicker in the PP group than in the normal group and the muscular layer of the lower part of the anterior uterine wall was very thin and adhered to the bladder ( Figure 1A). Ultrasonic examination showed penetrating placenta implantation, the disappearance of the posterior placental gap at the lower incision of the anterior uterine wall, and rich comb-like blood flow in the pernicious placenta previa group ( Figure 1B). Operative uterine and placental specimens showed that the placenta completely covered the internal os of the cervix was implanted in the myometrium and penetrated the serous layer ( Figure 1C). H&E staining results indicated that the myometrium was absent in the placental tissues of the percreta group, whereas the normal group had normal placental tissue ( Figure 1D). These results clarified the pathological features of the placental tissues for sequencing analysis.

| Different expression of coding genes, miRNAs and lncRNAs in PP
To determine the potential pathogenesis of PP, the expression profiles of genes (coding genes, miRNAs and lncRNAs) were determined by sequencing. A heatmap displays the differential expression of coding genes in the placental tissues of the placenta increta and normal groups (Figure 2A). The results showed that 469 coding genes were up-regulated in the PP group, whereas 321 coding genes were down-regulated, compared with the normal group (Table 2).
To confirm the accuracy of sequencing, qPCR was employed to determine the expression of four randomly selected coding genes (BTNL9, MAGEA4, ARHGEF28 and NR4A3). The expression trend determined by qPCR was consistent with the results of sequencing ( Figure 2B). Kyoto Encyclopedia of Genes and Genomes analysis indicated that up-regulated coding genes regulate the ErbB signalling pathway, bladder cancer and osteoclast differentiation ( Figure S1A), whereas down-regulated coding genes are involved in apoptosismultiple species, cytokine-cytokine receptor interaction and thiamine metabolism ( Figure S1B). miRNA sequencing showed that 178 miRNAs were up-regulated in the PP group, whereas 204 miRNAs were down-regulated, compared with the normal group ( Figure 2C; Table 2). The expression trends of miR-376c-3p, -655-3p, -3960 and -4492 were consistent between the sequencing and qPCR results ( Figure 2D). In total, 322 up-regulated lncRNAs and 219 down-regulated lncRNAs were identified in the PP group compared with the normal group ( Figure 2E; Table 2). The qPCR results of four randomly confirmed the accuracy of lncRNA sequencing ( Figure 2F). These results displayed the expression profile of lncRNAs, miRNAs and coding genes in PP.

| WGCNA of deregulated coding genes
To select the cluster of hub coding genes involved in the pathogenesis of PP, a WGCNA was carried out to divide the deregulated coding genes into several clusters based on a similar expression trend. The network of all gene clusters is displayed as a heatmap ( Figure 3A). The expression correlation was introduced as an important evaluation index to select functional modules. The ME black module and ME turquoise module were highly correlated with the pathogenesis of PP ( Figure 3B). Coding genes in the ME black were down-regulated and involved in cocaine addiction, fatty acid biosynthesis and melanogenesis ( Figure 3C,D). Coding genes in the ME turquoise were also down-regulated in PP and involved in protein processing in the endoplasmic reticulum (ER), lysosomes and N-glycan biosynthesis ( Figure 3E,F). These results indicated that coding genes in the ME black and ME turquoise play crucial roles in the pathogenesis of PP.

| Down-regulation of Wnt5A and MAPK13 in PP
WGCNA analysis suggested that both of ME black and ME turquoise play crucial roles in the pathogenesis of PP. Correlation analysis among coding genes was performed in ME black and ME  Figure 4E). Similar results were determined for MAPK13 expression in the normal and PP groups ( Figure 4F). Collectively, these results showed that down-regulated Wnt5A and MAPK13 in the ME turquoise are potential hub coding genes involved in the pathogenesis of PP.

| Down-regulation of lncRNA PTCHD1-AS and PAPPA-AS1 in PP
Next, we aimed to identify potential hub lncRNAs that correlated with Wnt5A and MAPK13. All deregulated lncRNAs in PP that correlated with Wnt5A and MAPK13 at the cut-off r value >0.5 were  Figure 7C).
These results were used to construct the lncRNA-miRNA-mRNA network.

| D ISCUSS I ON
In the present study, the results indicated that many lncRNAs, were collected as the normal group. According to the strict selection criteria, few samples were included in our study, which is the limitation of the present study.
In the present study, WGCNA was performed to select potential hub coding genes that were deregulated in the PP tissues.
The ME turquoise that contained down-regulated coding genes in PP tissues was used, and two potential hub coding genes, Wnt5A   Collectively, the present study demonstrated the crucial role of lncRNAs in the pathogenesis of PAS, which may provide a predictive biomarker and therapeutic target for PAS. However, further investigations are needed to expand the sample size and clarify the molecular mechanism involved.

ACK N OWLED G EM ENTS
Thanks to Chengdu Basebiotech Co., Ltd for providing assistance on bioinformatic analysis.

CO N FLI C T S O F I NTE R E S T
All authors declare that there are no conflicts of interest. Conceptualization (equal); funding acquisition (equal).

DATA AVA I L A B I L I T Y S TAT E M E N T
All data generated and/or analysed during this study are included in this published article.