Unveiling the immunomodulatory shift: Epithelial-mesenchymal transition Alters immune mechanisms of amniotic epithelial cells

Summary Epithelial-mesenchymal transition (EMT) changes cell phenotype by affecting immune properties of amniotic epithelial cells (AECs). The present study shows how the response to lipopolysaccharide of cells collected pre- (eAECs) and post-EMT (mAECs) induces changes in their transcriptomics profile. In fact, eAECs mainly upregulate genes involved in antigen-presenting response, whereas mAECs over-express soluble inflammatory mediator transcripts. Consistently, network analysis identifies CIITA and Nrf2 as main drivers of eAECs and mAECs immune response, respectively. As a consequence, the depletion of CIITA and Nrf2 impairs the ability of eAECs and mAECs to inhibit lymphocyte proliferation or macrophage-dependent IL-6 release, thus confirming their involvement in regulating immune response. Deciphering the mechanisms controlling the immune function of AECs pre- and post-EMT represents a step forward in understanding key physiological events wherein these cells are involved (pregnancy and labor). Moreover, controlling the immunomodulatory properties of eAECs and mAECs may be essential in developing potential strategies for regenerative medicine applications.


Gene-gene interaction network analysis reveals different controller genes
Using the data of KEGG pathways analysis, three gene-gene networks were built: LPS-eAECs network (network 1), LPS-mAECs (network 2) and a third network based on the LPS-mAECs vs. LPS-eAECs pairwise comparison (network 3).More in detail, network 1 consisted of 1,523 nodes and 7,880 edges with 14 connected components; network 2 consisted of 1,294 nodes, 11,792 edges, and 17 connected components; and network 3 consisted of 1,364 nodes and 7,625 edges grouped in 28 connected components.The topological parameters confirmed the scale-free nature of each network (Table 1).
Network 1 analysis identified 115 bottlenecks and 114 hubs (Table S6).Among them, kernel density estimation (KDE) analysis identified 4 local hubs: T cell receptor b chain V region (TRB), class II major histocompatibility complex transactivator (CIITA), T cell receptor beta variable 7-9 (TRBV7-9), and CD74 molecule (CD74), all genes involved in the adaptive immune response mechanism.Network 2 analysis identified 209 bottlenecks and 205 hubs (Table S6).Among the latter, 4 members of the signal transducer and activator of transcription (STATs) family are recognized as local hubs by KDE analysis: STAT6, STAT3, STAT5A, and STAT5B.Finally, Network 3 analysis identified 94 bottlenecks and 93 hubs (Table S6).Only one gene of those was identified as a local hub by KDE analysis, namely, the nuclear factor erythroid 2 like 2 (NFE2L2), which plays a pleiotropic role in the regulation of inflammation, autophagy, metabolism, and immune response mechanisms.Moreover, we analyzed the networks by means of the MCODE algorithm to identify different functional modules, which we further annotated in proper biological processes GO category by using ClueGO plugin.Network 1 was characterized by 9 functional modules (Figure 2) mainly associated with survival and immune stress response (fuchsia, pale blue, pink, green, and yellow modules) and antigen processing immune response processes (red and orange modules).
Modules related to RNA transport (purple) and metabolism (blue) were also identified (Table S7).Network 2 was characterized by 3 functional modules related to protein-mediated cell signaling response involving JAK-STAT and Notch signaling pathways (navy and dark purple modules).A module with genes playing a role in the regulation of synaptic vesicle processes (orange) was also identified (Table S7).Network 3 was characterized by 6 modules associated with cell defense (pink, pale blue, and blue modules).Modules related to infection response such as cardiac septum development, vesicle, and RNA trafficking processes (dark orange, dark red, and olive) were also identified (Table S7).

Candidate driver genes managing LPS response
To identify those genes that could be considered as potential ''drivers'' of the immune response induced by LPS response, we performed an integration of results from differential expression, topological network analysis, and module analysis (Figure 3; Table S8).
More in detail, in LPS-eAECs 34 genes displayed both hub and bottleneck node properties (Figure 3).In this category, RAN member RAN oncogene family (RAN), ras homolog family member A (RHOA), and eukaryotic translation initiation factor 2 subunit alpha (EIF2S1) also resulted to be differentially expressed (DEG) (Table S1).Two additional immune-related genes, T cell receptor b chain V region (TBR) and class II major histocompatibility complex transactivator (CIITA), also showed characteristics of local hubs.
The network analysis of the LPS-mAECs recognized 22 genes with hub-bottleneck properties.Among them, interleukin 6 (IL-6) and signal transducer and activator of transcription 2 (STAT2) also showed a DEG feature, whereas a role of local hub was also assigned to signal transducer and activator of transcription 5 (STAT5A) and signal transducer and activator of transcription 3 (STAT3) genes.
In addition, the comparative analysis between the two LPS-exposed cell populations networks highlighted 16 genes with hub-bottleneck features.Among them, STAT2 and epoxide hydrolase 1 (EPHX1) also displayed very high differential expression.Furthermore, KDE analysis identified within the hub-bottleneck category the nuclear factor erythroid 2 like 2 (NFE2L2) gene showing local hub features.These findings were summarized in Table 2.
Furthermore, the differential expression of some genes with specific topological features or exhibiting levels of fold change was validated by qRT-PCR (Figure S2).Briefly, in LPS-eAECs two local hubs (CIITA and CREB) were down-expressed while the B2M gene confirmed its overexpression.Instead, in the LPS-mAECs real-time qPCR confirmed the over-expression of the two bottleneck genes DDX58 and IL-6 (the latter also showed hub properties).Similarly, the NFLE2L gene that performed as hub in the comparison of the two LPS-stimulated populations was confirmed as over-expressed in the mAECs with respect to the eAECs (Figure S2).The table displays topological parameters computed, the node degree distribution, and correlation of node degree with clustering coefficient in the network.See also Table S6.
CIITA and Nrf2 regulate the response to LPS of AECs pre-and post-EMT In order to confirm whether CIITA and Nrf2 play key roles in the regulation of eAECs and mAECs immune response, respectively, transient gene depletion of these two factors by using specific siRNAs has been performed.The downregulation of either CIITA or Nrf2 was confirmed at mRNA (Figure 4A) and protein (Figure 4B) levels up to 48 h.Thereafter, siRNA-CIITA eAECs and siRNA-Nrf2 mAECs were exposed to LPS inflammatory stimulus for 1 h and serum-free conditioned media (CMs) were collected after 24 h.We investigated the immunomodulatory effect of CM derived from silenced eAECs and mAECs, respectively, through two different immune tests (Figures 4C-4F).Our data showed that eAECs and mAECs have two divergent immune behaviors; while the eAECs are mainly anti-inflammatory, the mAECs possess both anti-and pro-inflammatory properties.
In detail, the pro-inflammatory ability of the cells was determined by subjecting THP-1 macrophages to the collected CMs to assess the IL-6 production, used as a marker of their activation.In detail, CM derived from LPS-stimulated eAECs induced a significant increase of IL-6 production by THP-1 macrophages with respect to CM derived from non-stimulated eAECs (Figure 4C).Instead, the silencing of CIITA in eAECs significantly decreased IL-6 release with respect to non-silenced cells (Figure 4C).Interestingly, the transient depletion of CIITA in LPS-stimulated eAECs induced the highest IL-6 production, thus suggesting the anti-inflammatory role of this gene (Figure 4C).A similar trend of IL-6 production was recorded for CM derived from mAECs except for silenced mAECs stimulated with LPS (Figure 4D).In particular, THP-1 macrophages showed a significant increase of IL-6 production when subjected to CM derived from LPS-stimulated mAECs (Figure 4D).As with silenced eAECs, the transient depletion of Nrf2 in mAECs significantly decreased this release with respect to non-silenced cells (Figure 4D).Conversely, the CM derived from silenced mAECs did not impair the ability of mAECs to activate THP-1 macrophages when stimulated with LPS (Figure 4D).(C) LPS-mAECs vs. LPS-eAECs network, 6 modules: pink (glutathione metabolic process), pale blue (positive regulation of production of molecular mediator of immune response), blue (regulation of transcription from RNA polymerase II promoter in response to stress), dark orange (cardiac septum development), olive (translation initiation factor activity), dark red (synaptic vesicle exocytosis).See also Table S7.
Furthermore, the CMs collected from both eAECs and mAECs were added to phytohemagglutinin (PHA)-stimulated peripheral blood mononuclear cells (PBMCs) to evaluate CM effect on PBMCs proliferation.To this regard, the inhibition of PBMCs proliferation was used as an index of CM anti-inflammatory abilities.The obtained results demonstrated that CM derived from eAECs reduced PHA-stimulated PBMCs proliferation of 65% (Figure 4E).Moreover, CM derived from LPS-stimulated eAECs induced a significant decrease of PBMCs proliferation with respect to non-stimulated cells (Figure 4E).Of note, CIITA silencing dramatically impaired the eAECs-derived CM inhibitory effect on PBMCs proliferation with respect to non-silenced cells (Figure 4E).This trend was consistently maintained in the presence of LPS stimulation (Figure 4E).As with eAECs, CM derived from mAECs reduced PHA-stimulated PBMCs proliferation of 55% (Figure 4E).Besides, CM derived from LPS-stimulated mAECs induced a massive reduction of PBMCs proliferation with respect to non-stimulated cells (Figure 4E).Surprisingly, CM derived from silenced mAECs exhibited a similar inhibitory effect on PBMCs proliferation as recorded for CM collected from LPS-stimulated mAECs (Figure 4E).Finally, CM collected from silenced mAECs stimulated with LPS did not induce a comparable inhibition as recorded for silenced cells, although PBMCs proliferation was very low (Figure 4E).

DISCUSSION
Here, for the first time we reported that eAECs and mAECs possessed different immune behaviors and showed divergent immunomodulatory responses during inflammation, thus demonstrating that AECs phenotype may highly affect the immune response of these stem cells.In this context, we have recently shown that both eAECs and mAECs enhanced tendon regeneration in vivo by displaying distinct immunomodulatory profiles. 28In detail, the ability of these cells to convey the shift from pro-inflammatory to pro-regenerative responses in the host tissue differed in terms of macrophages polarization, cellularity, and blood vessel remodeling. 28Based on these premises, to further understand the molecular mechanisms underlying the cells' immune response, RNA-seq was used to compare the whole transcriptomics profiles in eAECs and mAECs.Subsequently, the adoption of network theory allowed identification of two different immune behaviors: eAECs mainly participate in the immune response by upregulating antigen-presenting circuit; on the other hand, mAECs mainly responded by upregulating cytokine and soluble mediators.Analysis of the gene involved in these two responses led to the identification of CIITA and Nrf2 as candidate driver genes, respectively, for eAECs and mAECs.Finally, the transient depletion of these genes in eAECs and mAECs demonstrated that they play a fundamental role for the respective immune responses: on the one hand, eAECs anti-inflammatory ability is dramatically compromised by the transient CIITA depletion.On the contrary, the contribution of mAECs to the immune system activation is highly impaired by the transient Nrf2 depletion, further suggesting that EMT may indeed represent a key event for the determination of AECs functions.
As epithelial cells, AECs play a fundamental role as a natural barrier between the organism and the environment (and fetus), developing, in turn, common strategies to regulate tolerance while preserving immunity against pathogens. 29Therefore, LPS stimulation of eAECs activates specific intracellular signaling belonging to different members of the TLR4/NF-kB/type I IFN pathway, identified as controllers of the interaction network.Of note, there is a group of genes endowed with antigen processing and presenting functions, such as CIITA, TRBV, rano class II histocompatibility antigen (a chain bL3-7 like and b chain like), b2-microglobulin (B2M), and RhoA.CIITA is a member of the  S8.NLR (nucleotide-binding and leucine-rich-repeat-containing) protein family and is recognized as the master regulator of MHC-II gene expression. 30Surprisingly, even though a master regulator of MHC-I gene has been identified (NLRC5), it was found that CIITA can indeed activate MHC-I expression, especially in cell lines like AECs 31,32 where those antigens are negative or have very low expression. 33,34Moreover, eAECs participate in the immune system activation through the upregulation of MHC-II proteins and all the complex of the proteins related to the antigen presenting process. 35,36To this regard, here we demonstrate that the transient depletion of CIITA dramatically impairs the anti-inflammatory abilities of eAEC, thus demonstrating its pivotal role in the immune response of these cells. 37onversely, the modulation of MHC genes seems not to be the primary mAECs mechanism in response to LPS.In fact, mAECs show the activation of a plethora of intracellular pathways specifically related to tyrosine kinase receptors (EGFR, PDGF, GHR, CSF3R, JAK/STATs) signaling and interleukins (IL-6, IL-12B, IL-10) production.In particular, the immune response of mAECs relies on the activation of STAT2 and IFN-stimulated genes, many of which encode proteins with antiviral, anti-proliferative, pro-apoptotic, or immunomodulatory functions. 38,39The upregulation of IL-6 may possibly be mediated by STAT2 itself in response to the classical NF-kB activators (e.g., LPS), as recently demonstrated. 40On the contrary, although STAT6 bears features of local hub, date hub, bottleneck, and DEG in mAECs, this gene is significantly downregulated.Interestingly, STAT6 is primarily involved in signaling activated by IL-4 and IL-13, which play roles in T helper 2 (Th2) differentiation, immunoglobulin class switching, macrophage activation, and MHC-II expression. 39This could partially explain the absent MHC-II gene modulation in mAECs.][43][44] However, the comparison between LPS-mAECs and LPS-eAECs points out a group of genes implicated in either EMT or immunomodulatory process (Figure 5).In particular, most of the driver genes belong to the Nrf2 pathway (e.g., EPHX1, GSTP1,and NFE2L2), a transcription factor that we previously described to be involved in EMT of AECs. 25,45,46Of note, mAECs immune response with respect to their epithelial counterpart (eAECs) is not strictly anti-inflammatory.To this regard, the dual role played by Nrf2 in leading either anti-inflammatory or inflammatory responses may potentially explain the results observed in mAECs. 47Indeed, it was extensively demonstrated that a functional cross talk exists between NF-kB and Nrf2 pathways and that can enhance or impair the inflammatory response, depending on the specific context and stimuli. 47Although from a functional perspective Nrf2 negatively controls the NF-kB pathway by multiple mechanisms, [48][49][50] our results indicated that LPS-mAECs preferentially activate both pathways.Moreover, the transient depletion of Nrf2 highly enhances their anti-inflammatory property and reduces the pro-inflammatory property, thus suggesting a potential involvement in the activation of inflammatory pathways, likely depending on the activation of the NF-kB pathway, even though further experiments are needed to demonstrate this hypothesis.
All together these results point out an intriguing scenario in which EMT plays a significant role in determining the immune response not only in isolated AECs but also in whole AM, where the two cell populations coexist.An interesting speculation is related to the possibility that AECs in AM could undergo EMT/MET to facilitate inflammation resolution by exploiting alternative strategies.At the same time, the two cell populations could cooperate either in damping inflammation or in remodeling AM during pregnancy and labor.In this scenario, epithelial or mesenchymal cells of AM would undergo EMT/MET upon LPS or other environmental stimuli (e.g., oxidative stress or pathogens).Notably, several authors reported that LPS can induce EMT via JAK/STAT and NF-kB signaling, [51][52][53][54][55] two of the most upregulated pathways in both eAECs and mAECs.However, further studies are necessary to explore this fascinating possibility and the relative consequences.involved in EMT in cancer [75][76][77] X X The main biological functions of each driver gene have been reported.The table also shows the topological and the expression features identified from computational analyses.See also Figure S2.

RNA extraction and library construction
RNA extraction was performed on each biological sample using RNeasy Mini Kit (#74106, QIAGEN) according to the manufacturer's instruction.RNA was quantified and checked for RNA integrity number (RIN) using the Agilent 2100 Bioanalyzer imposing a minimum value of 8.0 for processing.mRNA (poly-A) was isolated using oligo-dT beads (#20020595, Illumina) then purified using Dynabeads mRNA Purification Kit (#61006, Invitrogen).Libraries were prepared using the TruSeq Stranded mRNA Library Prep Kit (#20020594, Illumina) according to the manufacturer's instructions.

RNA sequencing and quality control
The sequencing reaction has been performed on an Illumina HiSeq 2500 platform in single-end for 75 cycles ensuring at least 30 million of reads per sample. 79To assess the quality of sequences, raw reads were subjected to quality check using fastQC tool (https://www.bioinformatics.babraham.ac.uk/projects/fastqc) and then trimmed by Trimmomatic software (https://bioinformaticshome.com/tools/rnaseq/descriptions/Trimmomatic.html#gsc.tab=0), 80imposing a minimum quality per base of 24 and a minimum length of reads equal to 36 bases (ILLUMINACLIP:TruSeq3-SE:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:24 MINLEN:36).Only those transcripts containing assembled regions were considered and only the longest alternative splicing transcripts from a specific gene were retained.The analysis was then conducted at gene level by assigning the corresponding gene to each transcript.

Gene ontology and KEGG pathways enrichment
For each gene, the annotation of GO terms for the three categories (http://geneontology.org) and the Kyoto Encyclopedia of Genes and Genomes (KEGG; https://www.genome.jp/kegg/)has been downloaded from the Ensembl database using the BioMart tool (https:// bioconductor.org/packages/biomaRt/). 83 DEGs were independently annotated in the three main GO ontologies (Biological process, Cellular component, and Molecular function) and successively it was carried out the enrichment analysis using the Fisher's test, implemented in the top GO R-package (https://bioconductor.org/packages/release/bioc/html/topGO.html) and Rgraphviz R-packages (https://www.bioconductor.org/packages/release/bioc/html/Rgraphviz.html).The enriched GO terms were also visualized by ReviGO (http://revigo.irb.hr/). 84Only those GO terms having an adjusted p value (q-value) % 0.01 were significantly enriched.The differentially expressed genes were also mapped to KEGG database using the clusterProfiler package (https://guangchuangyu.github.io/software/clusterProfiler/) 85to identify biological pathways enriched in each compared condition.The KEGG pathways enrichment analysis was carried out using a hypergeometric test implemented in the edgeR package (https://bioconductor.org/packages/release/bioc/html/edgeR.html). 86The pathways with a q-value %0.05 were recognized to be significantly enriched.

Gene act network building and analysis
The integration of differential expression with enriched KEGG pathways led to the building of a characteristic gene-act-network.Network was built using the KEGGgraph R-package (https://www.bioconductor.org/packages/release/bioc/html/KEGGgraph.html) 87; then to reduce complexity, the chemical compounds were removed from the final graph.For visualization, main component selection and topological features extraction Cytoscape 3.7.1 software was used. 88The hubs were identified as previously described, 89,90 by using the following equation: ND > m + s, where ND is the node degree, m is the mean node degree, and s is the node degree standard deviation.The identification of bottlenecks was carried out using the Cytoscape (Cytoscape Consortuim, http://www.cytoscapeconsortium.org/) plugin CytoHubba.It implements the following algorithm for bottleneck calculation: Ts is the shortest path tree rooted at node s.BN(v) = Ss˛V ps(v), where ps(v) = 1 if more than |V(Ts)|/4 paths from node s to other nodes in Ts meet at the vertex v; otherwise, ps(v) = 0 91 .

Modules extraction
The Cytoscape plugin Molecular Complex Detection (MCODE) 92 was used to extract and analyze clusters of highly interconnected nodes (modules) within the network.The significant modules with k-core >5 and nodes degree >5 were further selected for GO BP analysis to evaluate the biological processes where these modules involved in, which was performed using the Cytoscape ClueGO plugin.
Transient silencing of Nrf2 and CIITA with siRNA and conditioned media collection Transient gene depletion was performed by siRNAs using 50 nM ON-TARGET plus, small interfering RNA (siRNA) targeting Nrf2 (FE52L001000005, Dharmacon), CIITA (FE51L0010000050, Dharmacon) and Non-Targeting Control (hereafter referred as NTC) (FE5D0018101005, Dharmacon) and the TransIT-X2 Dynamic Delivery System (MIR-6000, Mirus) reagents, respectively, according to the manufacturers' recommendation.In detail, to improve the transfection, cells were grown until they reached 50-60% of confluency and then transient gene depletion was performed by using specific siRNAs.After 48 h of transfection cells were starved in serum free medium for 4 h and then they were immune activated with 1 mg/mL LPS for 1 h.Afterward, the medium was replaced with fresh serum free medium that was collected after 24 h (hereafter referred as conditioned medium; CM) for immune assays inclucing THP-1 macrophages IL-6 secretion and PBMCs proliferation assay.Cell pellets were collected to confirm transient gene depletion by Real-time qPCR and Western Blot.Target sequences of siRNAs have been reported in Table S9.

Real time qPCR
Total RNA was extracted with Direct-zol RNA kit (R2071, Zymo research) following the manufacturer's instructions.1 mg of total RNA was retrotranscribed using oligodT primers (BIO-38029, Bioline) and Tetro Reverse Transcriptase (BIO-65050, Bioline), following the manufacturer's instructions.qPCRs were carried out in triplicate using the SensiFAST SYBR Lo-ROX kit (BIO-94050, Bioline) on a 7500 Fast Real-Time PCR System (Thermo Fisher Scientific), according to the manufacturer's instructions.The following PCR conditions were used for all experiments: 95 C for 10 min, followed by 40 cycles at 95 C for 10 s and 60 C for 30 s. Relative quantification was performed by using the DDCt method.GAPDH (Glyceraldehyde 3-phosphate dehydrogenase) and YWHAZ (14-3-3 protein zeta/delta) were selected amongst housekeeping genes for gene quantification.Expression profiles were similar with both reference genes.Sequences of primers and conditions used in real-time qPCR are reported in Table S10.

Analysis of THP-1 macrophages IL-6 secretion
CM biological effect on human macrophages was performed by assessing the IL-6 release through an ELISA assay.In detail, THP-1-derived macrophages were treated with the CM diluted 1:1, with complete RPMI medium for 24 h.Thereafter, the supernatants were collected, centrifugated at 30003 g for 10 min at 4 C, filtered using a 0.2 mm Ministart sterile filter (#11740966, Sartorius), and stored at À80 C until usage.These cellular supernatants were used to quantify the concentrations of released IL-6 using the Human IL-6 Uncoated ELISA Kit assay (#88-7066-88, Invitrogen) according to the manufacturer's directions.The plates were read at 450 nm and the sensitivity of the used ELISA assay was in the range 2-200 pg/mL.

PBMCs proliferation assay
The immunomodulatory activities of eAECs and mAECs were analyzed on PBMCs by testing their proliferation.2 3 10^5 of PHA stimulated PBMCs were cultured for 72 h with CM derived from immune-activated eAECs amd mAECs (LPS stimulus) following targeting silencing (CIITA for eAECs; NRF2 for mAECs) and CM derived from non-tageting control (NTC for both eAECs and mAECs).PBMCs proliferation was assessed by using CellTiter96 Aqueous One Solution Cell Proliferation Assay following manifacturer's instruction (G3582, Promega).

Figure 1 .
Figure 1.KEGG pathway enrichment analysis (A-C) KEGG functional categories characteristics of LPS-eAECs (A), LPS-mAECs (B), and both LPS-stimulated cell populations in comparison (C).The Y axis represents the KEGG functional categories, while the X axis represents the number of DEGs.The KEGG subcategories are indicated at the end of each bar.See also supplementary information: Figure S1, Data S1, from Tables S1, S2, S3, S4, and S5.

Figure 2 .
Figure 2. The three gene-regulatory networks with their annotated modules identified by MCODE Cytoscape plugin and annotated by ClueGO (A) LPS-eAECs vs. eAECs network, 9 modules: red (positive regulation of myeloid leukocyte differentiation), orange (antigen processing and presentation of peptide antigen via MHC class I), yellow (structural constituent of cytoskeleton), green (malic enzyme activity), pale blue (cellular response to reactive oxygen species), blue (acetyl CoA biosynthetic process), fuchsia (positive regulation of NIK/NF-kappa B signaling), purple (translation initiation factor activity), pink (protein tyrosine kinase binding).(B) LPS-mAECs vs. mAECs network, 3 modules: navy (receptor signaling pathway via JAK/STAT), dark purple (regulation of Notch signaling pathway), orange (regulation of synaptic vesicles).(C)LPS-mAECs vs. LPS-eAECs network, 6 modules: pink (glutathione metabolic process), pale blue (positive regulation of production of molecular mediator of immune response), blue (regulation of transcription from RNA polymerase II promoter in response to stress), dark orange (cardiac septum development), olive (translation initiation factor activity), dark red (synaptic vesicle exocytosis).See also TableS7.

(
Continued on next page)

Table 1 .
Network topological parameters analysis for each experimental group

Table 2 .
Driver genes in LPS-stimulated AECs