DNA Microarray-Based Global Gene Expression Profiling in Human Amniotic Epithelial Cells Predicts the Potential of Microalgae-Derived Squalene for the Nervous System and Metabolic Health

In recent years, perinatal stem cells, such as human amniotic epithelial cells (hAECs), have attracted increasing interest as a novel tool of stem cell-based high-throughput drug screening. In the present study, we investigated the bioactivities of squalene (SQ) derived from ethanol extract (99.5%) of a microalgae Aurantiochytrium Sp. (EEA-SQ) in hAECs using whole-genome DNA microarray analysis. Tissue enrichment analysis showed that the brain was the most significantly enriched tissue by the differentially expressed genes (DEGs) between EEA-SQ-treated and control hAECs. Further gene set enrichment analysis and tissue-specific functional analysis revealed biological functions related to nervous system development, neurogenesis, and neurotransmitter modulation. Several adipose tissue-specific genes and functions were also enriched. Gene-disease association analysis showed nervous system-, metabolic-, and immune-related diseases were enriched. Altogether, our study suggests the potential health benefits of microalgae-derived SQ and we would further encourage investigation in EEA-SQ and its derivatives as potential therapeutics for nervous system- and metabolism-related diseases.


Introduction
In recent years, stem cells have been used extensively in high-throughput drug screening, from identifying target compounds to evaluating compounds at the preclinical stage [1][2][3]. In addition, the recent identification and isolation of human pluripotent adult stem cell lines from cord blood, adipose tissue, and bone marrow have offered opportunities in drug discovery and increased confidence in the mechanism of action of new targets [4]. However, although stem cell-based methods reduce the period of product development and rate of attrition of new therapeutics, their invasive extraction procedures, expensive cell reprogramming and maintenance procedures, and ethical constraints are the primary impediments to considering them as practical sources for drug screening [5].
In this regard, perinatal stem cells, such as human amniotic epithelial cells (hAECs), are gaining interest as a novel source of pluripotent stem cells (PSCs) with significant placed into 50 mL conical tubes with a pre-digestion buffer (CMF-HBSS, 0.5mM EGTA, Wako Pure Chemical Industries Ltd., Osaka, Japan, 20 mL) and were incubated at 37 • C for 10 min. Then the membrane was incubated in 0.05% trypsin-EDTA (Thermo Fisher Scientific, Waltham, MA, USA) at 37 • C for 40 min and placed on ice after incubation. Finally, Dulbecco's Modified Eagle Medium (DMEM, 10% FBS, 1% Penn-Strep, two volumes) was added and centrifuged at 200 g at 4 • C for 10 min to detach the epithelial cells from the membrane. The supernatant was discarded and the pellet was resuspended in 10 mL of DMEM. A 100 µm strainer (Merck EMD, Darmstadt, Germany) was used to remove small tissue aggregates, and clots from the cell suspension and the primary AECs were isolated. On 0.4% trypan blue dye (Dojindo Laboratories, Kumamoto, Japan) cell viability of AECs was found to be >90%. hAECs were maintained in the PromoCell Placental Epithelial Cell Basal Medium (PromoCell, C-26140).

hAEC Treatment
hAECs were treated with EEA-SQ at the concentration of 20 µg/mL for 7 days following the previous protocols of hAEC treatment [15,17,19,20]. The medium was changed with EEA-SQ every 48 h (total three times during the 7-day treatment period) following the manufacturer's protocol (ElplasiaTM, Kuraray Co., Ltd., Tokyo, Japan). Control samples were in the basal medium (PromoCell Placental Epithelial Cell Basal Medium), which was also changed every 48 h. Finally, RNA samples were collected for three biological replicates of each EEA-SQ-treated and control hAECs on day 7 and prepared for microarray analysis as reported in our previous studies [15,17,19,20]. We attempted time-dependent microarray analysis on hAEC using different natural compounds. We found that day 7 to day 10 were the best time points to study the early biological events induced by specific compounds. Before day 7, the effects of the compounds cannot be observed. After day 10, hAECs begin to determine their fate choice, and therefore, not suitable for functionality screening of the compounds.

RNA Extraction and Quantification
We used the Isogen kit (Nippon Gene Co., Ltd., Tokyo, Japan) to extract total RNA. RNA quantity and quality were determined using the NanoDrop 2000 spectrophotometer (ThermoScientific, Tokyo, Japan).

Microarray Gene Expression
We performed microarray analysis with Affymetrix's GeneAtlas ® System (Affymetrix Inc., Santa Clara, CA, USA). Firstly, total RNA samples were prepared from the three replicates of each control and EEA-SQ-treated hAEC sample (250 ng per RNA sample). Next, the amplified and biotin-labeled complementary RNA samples (cRNAs) were generated using the GeneChip 3 IVT PLUS Reagent Kit (Affymetrix Inc., Santa Clara, CA, USA). Then the cRNA samples were fragmented and labeled and then prepared for hybridization using the GeneAtlas ® Hybridization, Wash, and Stain Kit for 3' IVT Array Strips. Finally, the human genome array strips (HG-U219) were hybridized for 16 h at 45 • C on the GeneAtlas ® Hybridization Station, washed and stained with the GeneAtlas ® Fluidics Station, and scanned with the GeneAtlas ® Imaging Station.

Microarray Data Processing and Analysis
The raw image data were normalized following the robust multichip average (RMA) algorithm using the Expression Console Software (Affymetrix, Japan, URL: http://www. affymetrix.com, accessed on 30 November 2017). Subsequent gene expression analysis was carried out using freely available software Transcriptome Analysis Console (TAC) version 4 (Thermofisher Inc., Tokyo, Japan). Genes with a linear fold change >1.1 and p-value < 0.05 (one-way between-subjects ANOVA) were considered as differentially expressed genes (DEGs). The tissue specificity of the DEGs was determined using the Tissue-Specific Expression Analysis (TSEA) web tool (URL: http://genetics.wustl.edu/jdlab/tsea/, accessed on 14 October 2021). Tissue enrichment of the DEGs was identified by Fisher's exact test with the Benjamini-Hochberg correction and was expressed as the 'Specificity Index' thresholds (pSI) of varying stringency from 0.05 to 0.001. Genes with pSI < 0.0001 were the most specific genes enriched in a particular tissue [33,34].
HumanBase public database (URL: https://hb.flatironinstitute.org, accessed on 22 October 2021) and its analytical tool 'Module' was used to identify functionally clustered modules by the DEGs based on shared k-nearest-neighbors (SKNN) and the Louvain community-finding algorithm [36]. The SKNN-based strategy follows several steps to identify highly specific edges and to link genes that are highly likely to be part of the same functional cluster. Firstly, a user-selected tissue-specific network subset is created that contains the user-provided genes (nodes) and each edge between nodes, which are associated with a certain weight. Based on the original weight of the edges, a new weight is calculated that is equal to the number of k nearest neighbors. Then the top 5% of the edges, according to new edge weights, are selected to create the shared-nearest-neighbor tissue-specific network. Finally, the Louvain algorithm is applied to cluster this network into distinct functional modules. In order to stabilize the clustering, the Louvain algorithm is run 100 times and the score of cluster co-membership of each pair of genes/nodes equal to the fraction of times (out of 100) the pair was assigned to the same cluster is calculated. Genes with a co-membership score ≥0.9 are assigned to clusters. The resulting modules are tested for functional enrichment using genes annotated to Gene Ontology (GO) biological process (BP) terms and are presented with their resulting Q-value. The Q-value of each term within the modules is calculated using the one-sided Fisher's exact tests followed by Benjamini-Hochberg corrections. Similarly, tissue-specific functional modules were identified for both up and downregulated DEGs. Additionally, tissue and biological process-specific gene networks, tissue expression, and BP terms were retrieved for the top 19 up and downregulated DEGs.
Curated gene-disease association data were retrieved from the Comparative Toxicogenomics Database (CTD) (URL: http://ctdbase.org//, accessed on 22 October 2021) [37]. CTD contains two types of gene-disease associations-'curated' and 'inferred'. The curated associations are either extracted from the published literature or derived from the OMIM database using the NCBI Gene database. Inferred associations are based on CTD-curated chemical-gene interactions. We identified the diseases that are statistically enriched by the DEGs using the Set Analyzer tool that performs set-based enrichment analyses and pathway generation for user-provided genes. The significance of enrichment was calculated by the hypergeometric distribution adjusted by the Bonferroni method. MyGeneVenn tool was used to create Venn diagrams to compare our gene list to genes with direct relationships to specific disease terms (filtered for exact terms only). We checked for both inferred and curated gene-disease associations. Finally, genes with direct evidence of disease association were identified and presented according to their inference score. Inference scores are calculated for all inferred relationships and represent the similarity degree between the CTD chemical-gene-disease networks and a similar scale-free random network.
All data supporting the findings of this study are included within the article or its Supplementary Files. We have deposited the microarray data into the Gene Expression Omnibus (GEO) under the Accession Number: GSE188411.

Ethics Approval
The protocol was approved by the Ethical Review Committee of the University of Tsukuba Hospital (protocol approval code: H27-58, date: 10 July 2015). Informed written consent was obtained from the placenta donors.

Characteristics of Gene Expression Profiling in EEA-SQ-Treated hAECs
A total of 49,386 probe sets and 20,083 unique genes were identified on HG-U219 array strips. At day 7, 1841 genes were differentially expressed in the EEA-SQ treated hAECs compared to control hAEC (fold change > 1.1 and p-value < 0.05). Among them, 1256 DEGs were upregulated and 584 were downregulated. The scatter plot shows the signal intensities of the DEGs in EEA-SQ-treated and control hAECs (Figure 1a). Up and downregulated DEGs are shown in yellow and green dots, respectively. Bar graphs show the distribution of fold changes of up and downregulated DEGs ( Figure 1b).
As reported in our previous studies, we have also compared the day 7 EEA-SQtreated and control data sets with the day 0 control hAECs. Consistent with our previous findings [15,17,20], a comparison between day 0 and day 7 control hAECs (untreated cells) showed only early biological events related to cell differentiation, cell cycle arrest, and regulation of cell proliferation (data not presented). In contrast, a comparison between day 0 control and day 7 EEA-SQ-treated hAECs showed similar, but more significant, effects observed when we compared between day 7 EEA-SQ-treated and untreated (day 7) control hAECs. Therefore, in the present study, we restricted our analysis to only day 7 hAECs (treated vs. control) to avoid any exaggeration of the findings.

Tissue-Specific Enrichment Analysis of the Differentially Expressed Genes (DEGs)
Enrichment analysis of the expressions of the DEGs across 25 tissues was performed using the TSEA web tool (Figure 1c). Among the DEGs (n = 1841), 1703 genes were found to be expressed in at least one tissue. We found that the brain was the most significantly enriched tissue by the DEGs in all specificity index thresholds (pSI) varying from 0.005 to 0.0001. A total of 229 DEGs were brain-specific. Among them, 31 had the highest specificity threshold (pSI < 0.0001, p-value = 0.002), 35 had pSI < 0.001 (p-value = 0.001), 64 had pSI < 0.01 (p-value = 0.003), and 99 had pSI < 0.05 (p-value = 0.007). Table 1 shows the overlapped genes with pSI < 0.0001 and their functions. Brain-specific DEGs are involved in chemical synaptic transmission, axon development, neurotransmitter secretion, response to wounding, axon ensheathment, locomotion, negative regulation of neuron apoptotic process, astrocyte development, glutamate secretion, long-term synaptic potentiation, and positive regulation of neurogenesis.
The next significantly enriched tissue was the kidney. Several kidney-specific (pSI < 0.001) solute carrier genes were identified, responsible for ion transport. The skin and thyroid were also significantly enriched at pSI threshold < 0.0001. Skin-specific DEGs regulate several important biological processes, such as skin development, fatty acid metabolic process, keratinocyte differentiation, lipoxygenase pathway, linoleic acid metabolic process, regulation of water loss via skin. Thyroid-specific DEGs are involved in thyroxine biosynthesis, G protein-coupled receptors (GPCR) downstream signaling, and lipoprotein metabolism.
Adipose tissue and blood vessels were also significantly enriched; however, most of the overlapped genes were in a low specificity threshold (pSI < 0.05). Only two adipose tissue-specific genes had the highest stringency (pSI < 0.0001), which are involved in negative regulation of the lipid catabolic process, and insulin-stimulated glucose uptake in adipocytes. On the other hand, several DEGs were overlapped with testis-specific gene sets (pSI < 0.0001) but did not reach the statistical significance level (p-value = 1.0).

Gene Set Enrichment Analysis
Next, we identified the top enriched hallmark gene sets using the MsigDB (Figure 2a). Significance was measured as the false discovery rate (FDR q-value < 0.05; an analog of hypergeometric p-value following Benjamini-Hochberg correction. Hallmark gene sets 'summarize and represent specific, well-defined biological states or processes and display coherent expression'. There are 50 hallmark gene sets that condense information from over 4000 sets and, therefore, that have reduced noise and redundancy. We found that EEA-SQ significantly enriched several cell cycle-related genes sets, such as 'Genes involved in the G2/M checkpoint', 'Genes encoding cell cycle-related targets of E2F TFs' as well as 'Genes defining epithelial-mesenchymal transition'. EEA-SQ also significantly enriched gene sets of signaling pathways regulate the core cellular processes of proliferation, differentiation, and migration, such as 'Genes upregulated by activation of WNT/CTNNB1', 'Genes upregulated by activation of hedgehog signaling', and 'Genes upregulated by activation of Notch signaling'. Several inflammatory pathwayrelated gene sets were also regulated by EEA-SQ, such as 'Genes defining inflammatory response', 'Genes upregulated by PI3K/AKT/mTOR pathway activation', 'Genes upregulated by STAT5 in response to IL2 stimulation', 'Genes upregulated in response to TGFB1 , and 'Genes downregulated by KRAS activation'. . Significance was considered at false discovery rate (FDR) q-value < 0.05.

KEGG Pathways Enriched by EEA-SQ
Significantly enriched KEGG pathways (FDR q-value < 0.05) were identified on the MSigDB web tool (Figure 2b). Along with Wnt and Notch signaling pathways, EEA-SQ significantly enriched MAPK, calcium, chemokine, RIG-I-like receptor, neurotrophin, NOD-like receptor, ErbB, and Toll-like receptor signaling pathways. Interestingly, EEA-SQ also regulated adipocytokine and insulin signaling pathways. Other important pathways include long-term potentiation, melanoma, neuroactive ligand-receptor interaction, and vascular smooth muscle contraction.

Detection of Gene Clusters and Functional Module
We detected functionally clustered genes based on SKNN and the Louvain community-finding algorithm using the HumanBase public database, enabling network-based functional interpretation of the DEGs. Each module consists of tightly connected genes with a co-membership score ≥ 0.9. For each module, GOBP terms were detected for functional enrichment ( Figure 3).

KEGG Pathways Enriched by EEA-SQ
Significantly enriched KEGG pathways (FDR q-value < 0.05) were identified on the MSigDB web tool (Figure 2b). Along with Wnt and Notch signaling pathways, EEA-SQ significantly enriched MAPK, calcium, chemokine, RIG-I-like receptor, neurotrophin, NOD-like receptor, ErbB, and Toll-like receptor signaling pathways. Interestingly, EEA-SQ also regulated adipocytokine and insulin signaling pathways. Other important pathways include long-term potentiation, melanoma, neuroactive ligand-receptor interaction, and vascular smooth muscle contraction.

Detection of Gene Clusters and Functional Module
We detected functionally clustered genes based on SKNN and the Louvain communityfinding algorithm using the HumanBase public database, enabling network-based functional interpretation of the DEGs. Each module consists of tightly connected genes with a co-membership score ≥ 0.9. For each module, GOBP terms were detected for functional enrichment ( Figure 3). . DEGs in M2 also enriched BP terms related to the development and morphogenesis of the renal tubule, urogenital system, blood vessels, and cardiovascular system (Supplementary File S1). DEGs in M3 (n = 427) significantly enriched regulation of body fluid levels (GO:0050878), long-chain fatty acid metabolic process (GO:0001676), water homeostasis (GO:0030104), and regulation of blood circulation (GO:1903522). Global functional modules and tissue-specific functional modules, and enriched BP terms in each module, have been described in Supplementary File S1.

Tissue-Specific Biological Functions
We performed tissue-specific network-based functional interpretation of the upregulated ( Figure 4) and downregulated ( Figure 5) DEGs using the HumanBase web tool. We identified functional modules specific to adipose tissue, blood vessels, brain, colon, heart, kidney, liver, lung, muscle, skin, and thyroid, considering significantly enriched tissues by the DEGs (Figure 1c). . DEGs in M2 also enriched BP terms related to the development and morphogenesis of the renal tubule, urogenital system, blood vessels, and cardiovascular system (Supplementary File S1). DEGs in M3 (n = 427) significantly enriched regulation of body fluid levels (GO:0050878), long-chain fatty acid metabolic process (GO:0001676), water homeostasis (GO:0030104), and regulation of blood circulation (GO:1903522). Global functional modules and tissue-specific functional modules, and enriched BP terms in each module, have been described in Supplementary File S1.

Tissue-Specific Biological Functions
We performed tissue-specific network-based functional interpretation of the upregulated ( Figure 4) and downregulated ( Figure 5) DEGs using the HumanBase web tool. We identified functional modules specific to adipose tissue, blood vessels, brain, colon, heart, kidney, liver, lung, muscle, skin, and thyroid, considering significantly enriched tissues by the DEGs (Figure 1c). We found that upregulated DEGs showed mostly brain-specific functional enrichment ( Figure 4). Except for the inflammatory pathways (negative regulation interleukin-6 production (GO:0032715), lymphocyte-mediated immunity (GO:0002449), and acute inflammatory response (GO:0002526)), almost all the global BP terms were also brain-specific. Other important brain-specific BPs were positive regulation of neuron differentiation (GO:0045666), positive regulation of neurogenesis (GO:0050769), positive regulation of neuron projection development (GO:0010976), and positive regulation of nervous system development (GO:0051962).
There were also adipose tissue, skin and thyroid-specific functional enrichment by the upregulated DEGs. Adipose tissue-specific top BP terms were the transmembrane receptor protein tyrosine kinase signaling pathway (  We found that upregulated DEGs showed mostly brain-specific functional enrichment ( Figure 4). Except for the inflammatory pathways (negative regulation interleukin-6 production (GO:0032715), lymphocyte-mediated immunity (GO:0002449), and acute inflammatory response (GO:0002526)), almost all the global BP terms were also brain-specific. Other important brain-specific BPs were positive regulation of neuron differentiation (
Interestingly, adipose tissue-specific interactions were broadly similar for both top downregulated and top upregulated DEGs. The top downregulated DEGs were functionally interconnected to (adipose tissue-specific) CAV1, CD36, IL6, ADIPOQ, and AGT. Other adipose tissue-specific functionally connected genes were serglycin (SRGN), selenoprotein P (SELENOP), adaptor related protein complex 1 subunit sigma 2 (AP1S2), nuclear receptor coactivator 7 (NCOA7), annexin A1 (ANXA1), TIMP metallopeptidase inhibitor 1 (TIMP1), and CD44 molecule (CD44). The top enriched functions represented within the adipose tissue-specific network were response to wounding, negative regulation of cytokine secretion, low-density lipoprotein particle clearance, ERK1 and ERK2 cascade, and intrinsic apoptotic signaling pathway. The lists of top up and downregulated DEGs and their functions are presented in Tables 2 and 3.    Tricarboxylic acid metabolic process, drug metabolic process, energy derivation by oxidation of organic compounds, cellular respiration

Discussion
In the present study, we have performed integrated transcriptome analysis of EEA-SQ-treated hAECs. Our findings suggest the prospects of EEA-SQ in nervous system diseases and adipose tissue dysfunctions.
Natural resources are chemically diverse and are able to modulate several biological functions simultaneously. DNA microarrays have become a necessary tool to investigate the mechanisms of action and drug development of natural bioactive compounds [38,39]. While the conventional gene expression profiling are optimized for single gene analysis, microarrays generate gene expression data on a genome-wide scale and therefore permit systematic approaches to the biological discovery of regulatory mechanisms and biochemical pathways [40]. However, even employing high throughput tools like DNA microarrays may not improve the success rate for drug discovery and development as the essential biological elements of drug-target interactions are mostly lost in the monolayer Next, we identified the genes in our gene list that have inferred or curated association to nervous system diseases, namely neurodegenerative diseases, premature aging, depression and anxiety, and adipose tissue dysfunctions, such as obesity and diabetes mellitus (DM). Curated gene-disease associations are determined by CTD and OMIM curation and represent a more robust association compared to the inferred gene-disease associations, which are obtained via CTD-curated chemical-gene interactions. We selected the genes that have only direct relationships to the disease terms and filtered out the hierarchical relationship genes. We found that 60 DEGs from our gene list were directly associated with neurodegenerative diseases, 9 DEGs to depressive orders, 1 DEG to anxiety disorders, 14 DEGs to obesity, and 37 DEGs to DM (Figure 8b).
From curated gene-disease associations, we identified several genes that may be responsible for the mechanism of action of EEA-SQ. Next, we examined the inference score of those genes to select the DEGs with specific disease associations. The inference score represents the similarity degree between CTD chemical-gene-disease networks and scale-free random networks. The higher the score, the more likely the inference network has atypical connectivity. Brain-derived neurotrophic factor (BDNF), fibroblast growth factor 2 (FGF2), fibroblast growth factor receptor 1 (FGFR1), paraoxonase 1 (PON1), neuropeptide S (NPS), glial fibrillary acidic protein (GFAP), calcium/calmodulin-dependent protein kinase kinase 1, alpha (CAMKK1), nucleotide-binding oligomerization domain containing 1 (NOD1), serine hydroxymethyltransferase 1 (soluble) (SHMT1), and perilipin 1 (PLIN1) were the most important DEGs that have a curated association to a disease (either therapeutic target or marker of a disease) or inferred association with a very high inference score (Figure 8c). Their expression profile in EEA-SQ-treated hAECs, BPs, and disease association are described in Supplementary File S3.

Discussion
In the present study, we have performed integrated transcriptome analysis of EEA-SQtreated hAECs. Our findings suggest the prospects of EEA-SQ in nervous system diseases and adipose tissue dysfunctions.
Natural resources are chemically diverse and are able to modulate several biological functions simultaneously. DNA microarrays have become a necessary tool to investigate the mechanisms of action and drug development of natural bioactive compounds [38,39]. While the conventional gene expression profiling are optimized for single gene analysis, microarrays generate gene expression data on a genome-wide scale and therefore permit systematic approaches to the biological discovery of regulatory mechanisms and biochemical pathways [40]. However, even employing high throughput tools like DNA microarrays may not improve the success rate for drug discovery and development as the essential biological elements of drug-target interactions are mostly lost in the monolayer cell cultures [41]. In the present study, we used 3D hAEC spheroid to understand better the effect of EEA-SQ in a physiologically relevant context.
Previously we showed that different natural compounds could direct the differentiation of hAECs towards multiple cell lineages. A caffeoylquinic acid derivative-TCQA [20] and a caffeic acid ester-rosmarinic acid [15] enhanced neural cell differentiation in hAECs, whereas a flavonol, isorhamnetin, induced hepatic-lineage specific [21] and an anthocyanin CY3G induced adipocyte differentiation [19] in hAECs. We noticed that the bioactivities of natural compounds could broadly be predicted from the enriched cell types by differentially expressed genes. In the present study, we found that the DEGs by EEA-SQ treatment in hAECs could significantly enrich brain tissue, i.e., EEA-SQ primarily regulated brainspecific gene expression. Additionally, adipose tissue, blood vessel, colon, kidney, skin, thyroid tissues were enriched (Figure 1c). However, when we performed tissue-specific functional enrichment analysis, only brain, adipose tissue, and skin-specific biologic functions (BPs) remained relevant (Figures 4 and 5). The hallmark gene sets identified cell cycle-related functions, and core pathways involved in cell proliferation and differentiation, such as Wnt, Notch, and hedgehog signaling ( Figure 2a); however, KEGG pathway analysis could additionally identify the MAPK pathway, chemokine, calcium and neurotrophin signaling pathways, long-term potentiation and neuroactive ligand-receptor interaction (Figure 2b). When we performed tissue-specific functional network analysis, it became evident that EEA-SQ mainly significantly regulated nervous system development, neuron differentiation, neurogenesis, and neuron projection development. Brain-specific DEGs are also involved in neurotransmitter secretion and antiapoptotic and anti-inflammatory functions in the brain (Table 1). When we limited our analysis to only the top regulated DEGs (highest fold change), we found that the top upregulated DEGs are primarily expressed in the brain and are functionally related to multiple brain-specific and neuron-differentiationspecific genes (Figure 7). For example, the top upregulated DEGs are functionally associated with BRCA1, which is highly expressed in adult neurogenic areas and regulates distinct apoptotic functions in neural progenitors [42]. Further gene-disease association analysis could identify some target genes of EEA-SQ for its therapeutic potential in nervous system-related diseases, such as BDNF, FGF2, GFAP, PON1, and NPS. BDNF is a neurotrophic factor that promotes neuronal survival, growth, and maturation [43][44][45]. FGF2, a multi-functional growth factor, is expressed in the subventricular zone (SVZ) and the subgranular zone (SGZ) of the hippocampal dentate gyrus and regulates neural stem and progenitor cells and, therefore, has been implicated in adult neurogenesis [46]. GFAP is expressed in the astrocytes and is involved in cell communication and the functioning of the blood-brain barrier in the central nervous system. On the other hand, PON1 is associated with neurobehavioral development [47,48]. A neuropeptide NPS increases the excitatory transmitter glutamate secretion and produces anxiolytic-like effects by reducing fear responses [49]. All these findings align with our previous studies on EEA-SQ that reported antidepressant-like, neuroprotective, and cognition-enhancing effects of EEA-SQ in vivo [31,32].
Additionally, some critical adipose tissue-specific functions were observed in EEA-SQtreated hAECs. Significantly enriched KEGG pathways include adipocytokine and insulin signaling pathways and tissue-specific functional analysis reveals positive regulation of MAPK cascade. Furthermore, functional network analysis of the top DEGs identified several functionally connected genes in adipose tissue, including PPARG, SERPINE1, CD36, CAV1, CAV2, and ADIPOQ. PPARG is the key regulator of adipogenesis, is involved in fatty acid metabolism and triglyceride storage, and suppresses the production of inflammatory mediators in macrophages [50]. SERPINE1 is a unique insulin-sensitizing adipocytokine and, like PPARG, is implicated in a wide range of metabolic diseases, including DM, obesity, and fatty liver diseases [51]. CD36 is a marker of human adipocyte progenitors [52], CAV1 and CAV2 are important regulators of adipose tissue homeostasis [53], and ADIPOQ is observed exclusively in mature fat cells [54]. Indeed, the top enriched BPs represented in the adipose tissue-specific network included antiinflammatory, lipoprotein particle receptor biosynthetic process, and lipoprotein transport as well as lipid droplet formation. Lipid droplets in the form of foam cells derived from macrophages are a hallmark for the progression of atherosclerosis [55], which explains the enrichment of cardiovascular diseases by the DEGs. We also found that EEA-SQ significantly upregulated the expression of PLIN1, a regulator of lipid storage and lipolysis [56]. It is worth noting that in our previous study we reported that SQ and its derivative improved the metabolism of adipocytes differentiated from diabetic adipose-derived stem cells by enhancing energy homeostasis and insulin sensitivity and also prevented excessive lipogenesis [57].
Previously, we have reported anti-inflammatory effects of Aurantiochytrium strains in murine macrophages [29,30]. In hAECs, EEA-SQ also significantly enriched antiinflammatory hallmark gene sets and KEGG pathways. However, functional analysis could not detect any tissue-specific anti-inflammatory function. Therefore, we assume that EEA-SQ may regulate the systemic anti-inflammatory function instead of the tissue-specific mechanism of inflammation resolution.
Additionally, our TSEA analysis showed that skin is one of the significantly enriched tissues (Figure 1c). Functional module analysis (Figure 3) also revealed enrichment of the dermatan sulfate proteoglycan biosynthetic process (GO:0050651) and exosomal secretion (GO:1990182). It has been reported that alteration of dermatan sulfate proteoglycan biosynthesis in the fibroblasts of skin results in severe skin structure disruption [58,59]. Exosomes are made of a lipid bilayer structure containing lipids, nucleic acids and proteins and are secreted by most cell types. Previous studies reported the efficacy of exosomes in a number of skin defects, including skin aging, skin inflammation, atopic dermatitis, and pigmentation [60][61][62][63]. Also, tissue-specific functional analysis showed that the long-chain fatty acid metabolic process, water homeostasis, and antiinflammatory functions were significantly enriched in skin (Figure 4). SQ is very popular in cosmetics products as it maintains skin's moisture barrier and hydration [25]. Our findings suggest that EEA-SQ may also have similar beneficial properties for skin.
Although our study explores multiple beneficial prospects of EEA-SQ, our results at present are observational in nature at the whole-genome transcriptome level. However, some of our observations are consistent with our previous animal-model studies [31,32] that reported antidepressant and neuroprotective effects of EEA-SQ. Further in vitro study focusing on the time-and dose-dependent effects of EEA-SQ on hAEC morphology and phenotypes are required to confirm the potential functionalities of EEA-SQ observed in the present study.

Conclusions
Altogether, our findings suggest the potential health benefits of EEA-SQ, notably in neurological and neurodegenerative diseases and metabolic disorders with adipose tissue dysfunction. Thus, our work may promote interest in and encourage investigation of EEA-SQ and its derivatives as potential therapeutics for nervous system-and metabolismrelated diseases.

Institutional Review Board Statement:
The study was conducted according to the guidelines of the Declaration of Helsinki, and the protocol was reviewed and approved by the Ethical Review Committee of the UNIVERSITY OF TSUKUBA (approval code: H27-58, approval date: 10 July 2015).

Informed Consent Statement:
Informed consent was obtained from all the mothers who donated the placenta. Data Availability Statement: All data generated or analyzed during this study are included in this published article and its Supplementary Information Files. Microarray data are deposited in the Gene Expression Omnibus (GEO) under Accession Number: GSE188411 (https://www.ncbi.nlm.nih.gov/ geo/query/acc.cgi?acc=GSE188411/, accessed on 13 November 2021).