Comprehensive analysis of differential immunocyte infiltration and the potential ceRNA networks during epicardial adipose tissue development in congenital heart disease

Background To detect the development, function and therapeutic potential of epicardial adipose tissue (EAT); analyze a related gene expression dataset, including data from neonates, infants, and children with congenital heart disease (CHD); compare the data to identify the codifferentially expressed (DE) mRNAs and lncRNAs and the corresponding miRNAs; generate a potential competitive endogenous RNA (ceRNA) network; and assess the involvement of immunocyte infiltration in the development of the EAT. Methods Multiple algorithms for linear models for microarray data algorithms (LIMMA), CIBERSORT, gene-set enrichment analysis (GSEA), and gene set variation analysis (GSVA) were used. The miRcode, miRDB, miRTarBase, and TargetScan database were used to construct the ceRNA network. The Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses of the DE mRNAs were performed. Results Thirteen co-DE mRNAs and 47 co-DE lncRNAs were subsequently identified. The related categories included negative regulation of myoblast differentiation, regulation of ion transmembrane transport, and heart development, which were primarily identified for further pathway enrichment analysis. Additionally, the hub ceRNA network in EAT development involving MIR210HG, hsa-miR-449c-5p, and CACNA2D4 was generated and shown to target monocyte infiltration. Conclusion These findings suggest that the pathways of myoblast differentiation and ion transmembrane transport may be potential hub pathways involved in EAT development in CHD patients. In addition, the network includes monocytes, MIR210HG, and CACNA2D4, which were shown to target the RIG-I-like receptor signaling pathway and PPAR signaling pathway, indicating that these factors may be novel regulators and therapeutic targets in EAT development.

in the myocardium and vascular arteries [1]. There are key metabolic differences between the EAT and other adipose tissue depots. Importantly, a systematic review from Iacobellis's group revealed that the rate of fatty acid synthesis, release, and breakdown was significantly higher in EAT than other visceral fats, and EAT exhibited lower glucose utilization as well [2]. Additionally, EAT appears to be rich in saturated fatty acids and metabolically related proteins and may have protective effects on the myocardium and vascular arteries by supplying the energy substrates for use in the ischemic myocardium and against elevated levels of FFAs in the myocardium and vascular damage to the arteries [3]. The higher lipolytic potential of EAT may also help meet the increasing myocardial energy demands, especially in the case of exercise-related metabolism or ischemic conditions [1,4,5]. Although fatty acids may contribute to the major energy substrates that are used to maintain cardiac metabolism and contractile function, the uptake of FFAs also leads to the impairment of mitochondrial structure and function as well as accumulation of toxic metabolites, including membrane potential disorder, Ca 2+ overload, and cytochrome c (Cyt-C) release; these changes induce oxidative stress and disrupt metabolic pathways [6]. Anatomically, EAT shares a common blood supply with the underlying myocardium without the fascia-like structure [5]. Researchers increasingly believe that EAT contributes to pathophysiological processes related to coronary atherosclerosis, diastolic dysfunction, atrial fibrillation, and progressive ventricular remodeling, exerting exocrine and paracrine effects on the heart vasculature [4,7,8]. Recently, data have increasingly showed that EAT contains a high density of lymphoid clusters, and it is generally accepted to be a source of immunomodulatory and inflammatory mediators [9]. Of note, there is a mutually reinforcing relationship in the pathological processes of EAT formation and cardiac remodeling.
For adult patients, Patel et al. revealed that EAT could secrete multiple adipokines, including apelin, adiponectin, adrenomedullin, angiotensinogen, leptin, and visfatin, coupled with the potential paracrine function of regulating inflammation, immune response, coronary plaques, ventricular remodeling, and vascular smooth muscle cell growth and migration [5]. Sakamoto et al. discovered a significant association between immuno-inflammatory mediators and coronary artery calcification, which may reflect immuno-inflammation mechanisms related to the development of cardiovascular remodeling [7]. Interestingly, inflammatory cytokines, inflammatory adipokines, and M1/M2 macrophages significantly accumulate in EAT, aggravating ischemic conditions as well as the damage due to myocardial infarction [7]. Park et al. analyzed the evidence from clinical studies and clearly demonstrated that EAT thickness could predict the occurrence of metabolic syndrome and coronary atherosclerosis in an Asian population and noted that it was stronger in patients with BMI < 27 kg/m 2 [10]. Additionally, for pediatric patients, the nutrition demand and morphology of organ development in the infant stage is critical for life. Ojha et al. found that there are unique gene expression patterns in EAT tissues in a comparison of thermogenic gene transcripts of in various developmental stages (~ 10%) [11]. Chechi et al. found that the expression of thermogenic genes significantly correlated with oxidative stress-and immune-related pathways, also presenting an important role in microenvironmental regulation within EAT tissue and local crosstalk between the cardiomyocytes and interstitial cells [12,13].
In this scenario, EAT, an immune and endocrine organ, may function as a key inflammatory mediator for cardiovascular system development. Here, the insulin-like growth factor 1 receptor (IGF1R) was shown to be stimulated after myocardial damage and then promotes the differentiation of epicardial cells into adipocytes to protect cardiomyocytes and ensure the survival of cardiac progenitors [14]. In contrast, EAT serves as a secondary lymphoid organ that regulates the innate and adaptive immune responses and thus locally and systemically promotes immune cell proliferation, activates inflammatory receptors, and recruits lipid deposition across the vessel wall [9]. In this study, we performed a systematic microarray analysis to reveal the gene expression changes, inflammatory cell infiltrates, and noncoding RNA-immunocyte coregulation networks that may help us identify novel biomarkers and therapeutic targets.

Data sources and processes
We performed a differential expression (DE) analysis based on "affy" package of R software that was curated by data processing principles includes expression data background correction, log2 transformation, data normalization and linear model and DE analyses. Here, the data normalization, missing value procession, and powerful DE analyses were performed by a function of the robust multi-array average (RMA), K-nearest neighbor (KNN), and linear models for microarray data algorithms (LIMMA), respectively [15]. We defined the DE threshold adjusted p-value, corrected by the Benjamini-Hochberg method, as less than 0.05 and |log 2 fold-change (FC)| greater than 1.5. Herein, these genes were identified as having significantly increased (log 2 FC > 1.5) or decreased (log 2 FC < − 1.5) expression. Figure 1 illustrates the schema of the analysis flow, and the statistical analysis was performed with the R (version 3.6.1) and Perl (version 5.30.1) toolkit.

Detection of immune infiltration subtypes
In this routine, we applied the CIBERSORT algorithm (https ://ciber sort.stanf ord.edu/) to estimate the burden of immune infiltration subtypes in the EAT development process, based on a linear support vector regression (SVR) algorithm with gene expression data with systems-level gene expression data in the R toolkit [16].
We chose the hub cell-subtypes with the parameters that were set as follows: (1) all the gene expression data were normalized by RMA algorithm; (2) the perm parameter was 1000; (3) the statistical differences between groups had p-value < 0.05; (4) multiple comparisons among the neonatal, infant and child groups were tested by one-way ANOVA, and the p-value was adjusted by the Bonferroni method.

Gene set bioinformatics enrichment analysis
Gene Ontology (GO) term enrichment analyses were carried out separately for the DE genes using the MetaScape database (http://metas cape.org/), which is a webbased analytics tool used for gene functional enrichment analysis [17]. Of note, the major GO terms of the DE genes in biological processes, molecular functions, and cellular components were illustrated. The parameter settings were as follows: (1) min overlap is 3; (2) the cutoff p-value = 0.05; (3) the min enrichment score = 1.5.

Immunocyte infiltration analysis
Differential immunocyte Step 1 Step 2 Step 3 Multiple comparisons Based on the time series-related Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway and biological process term enrichment, a gene-set enrichment analysis (GSEA; http://softw are.broad insti tute.org/gsea/index .jsp) was performed [18]. The parameters were as follows: (1) the threshold for significant enrichment was p-value < 0.05; (2) the permutation test number = 1000; (3) the gene profile was defined as continuous phenotypes with time series experiments for EAT development and defined as class labels in response to the expression level of hub mRNAs/immunocyte infiltration; (4) enrichment statistic was used for weighting; (5) the metric for ranking genes selected "Pearson" for time series analysis and defined the "Signal2Noise" for categorical analysis.

Identification of the downstream regulatory pathways of target genes
After grouping by the median expression value of the target gene, lncRNA, and infiltrating immunocytes, we applied the gene set enrichment analysis (GSEA) to detect the downstream pathways. The normalized enrichment score (NES) was obtained via cluster analysis based on the Pearson statistical method [18]. In addition, after identification of the potential regulatory pathways, the gene set variation analysis (GSVA) package, based on a nonparametric unsupervised method, was applied to explore the expression of hub genes in each sample in relation to the hub ceRNA network and candidate immunocytes in our study [24].

Validation of candidate immunocytes, regulatory networks and hub genes
Finally, we also verified the landscape of immunocyte infiltration, regulatory networks and hub gene expression related to EAT development based on the GSE115799 dataset for Ovis aries. The CIBERSORT algorithm was used to compare the EAT sample immune infiltration between 1 and 4 weeks of age. Furthermore, the candidate pathway variation was assessed by the GSVA package. The mRNA and lncRNA expression values between the two groups were compared with powerful DE analyses with the LIMMA algorithm.
Of note, 46 EAT samples were obtained from cardiopulmonary bypass cardiac surgery of pediatric patients with congenital heart disease (CHD). Here, for the comparison among the neonate, infant, and child EAT samples, to illustrate the transcription and expression profiles in response to EAT development, we performed a powerful differential analysis across the comparison groups. The results showing that 270, 889, and 1069 DE mRNAs were detected in the comparisons of infant vs child, infant vs neonate, and neonate vs child, respectively  . 1a) and biological process terms of Gene Ontology (GO) (Fig. 1b)  lncRNAs were identified to construct the heatmap. A cluster analysis based on the Pearson method presented a significant intergroup difference (Fig. 2c, d).

Immunocyte subtype detection
As shown in Fig. 3d, Table S4). As a result, the immunocytes of monocytes, as a hub immune cell fraction, were quantified from the gene expression matrix derived from the EAT samples (Fig. 3e).

CeRNA network construction
Our findings, including (1) 1742 links, were calculated based on the 47 co-DE lncRNAs and targeted miRNAs, (2) and 27 lncRNAs and 7 miRNAs were selected after removing the duplicate nodes, (3) and 4 targeted mRNAs were selected after determination of the miRNAs targeting mRNAs and co-DE mRNAs in multiple group comparisons. As a consequence, the ceRNA network including 27 lncRNAs, 7 miRNAs, and 4 mRNAs was accordingly constructed for EAT development for CHD patients ( Fig. 4a; Additional file 5: Table S5). A significant coefficient of correlation linking the infiltrating value of monocytes and the expression value of CACNA2D4 was identified (Pearson's correlation coefficient = 0.62, p-value = 4.32E−06) (Fig. 4b). Subsequently, a subnetwork including CACNA2D4, hsa-miR-449c-5p, C7orf13, PTEN, LINC00478, WDFY3.AS2, WHSC1L1, and MIR210HG was identified as a core regulatory network (Fig. 4d). The cluster analysis results of these biomarkers showed significant discrimination in response to different developmental stages of EAT (Fig. 4c). Then, to further elucidate the lncRNA-miRNA-mRNA regulatory mechanism in response to monocyte infiltration, we calculated Pearson's correlation coefficient values for these biomarkers. A high-density correlation module including monocytes, MIR210HG, and CACNA2D4 was generated ( Fig. 4e; Additional file 6: Table S6).

GSEA and GSVA analysis
The GSEA analysis was applied to elucidate the potential regulatory pathways after grouping by the median value.
The cluster analysis of the pathway NES value showed that the PPAR signaling pathway was significantly correlated with higher monocyte infiltration, higher CAC-NA2D4 expression, and higher MIR210HG expression ( Fig. 5a and Additional file 7: Table S7). In contrast, the RIG-I-like receptor signaling pathway was primarily associated with lower monocyte infiltration, lower CAC-NA2D4 expression, and lower MIR210HG expression (Fig. 5b).

Verification of immunocytes, regulatory maps and hub genes in an Ovis aries model
After analysis, the immuno-infiltration landscape is presented in Fig. 6a (Additional file 4: Table S4). Additionally, compared with that of the 7 days of age group, the infiltration level of monocytes was decreased in the 28 days of age group (Fig. 6b). The EAT transcriptomic levels of CACNA2D4, C7orf13, PTEN, LINC00478, WDFY3-AS2, TMEM72-AS1, WHSC1L1, and MIR210HG were compared between the 7 and 28 days of age groups, and the results showed that both the CACNA2D4 and MIR210HG genes were upregulated in the 28 days of the age group compared with the 7 days of age group (Fig. 5c,  d). Based on GSVA analysis, the PPAR signaling pathway (variation score p-value = 0.035) and the RIG-I like receptor signaling pathway (variation score p-value = 0.046) were verified and are thus considered important regulatory pathways/networks in EAT development (Fig. 6e).

Discussion
At present, EAT is believed to be a large secretosome that has major roles in regulating the physiological and pathophysiological processes in cardiovascular disease. Furthermore, EAT exerts a potent effect on the paracrine regulation of vascular atherosclerosis susceptibility, myocardial fibrosis, and arrhythmia based on functional and anatomic proximity. In sum, EAT has been associated with thermogenesis, immune-inflammation regulation, cardiometabolic activity, and adrenergic receptor distribution and thus has been a research hot spot in recent years. For the first time, we integrated human and O. aries microarray-based profiles of expression data to determine the immunocyte infiltration-and the candidate immunocyte-correlated regulatory ceRNA network during the development of EAT through early life. Further extensions of the bioinformatics tools presented herein were used to identify the molecular mechanisms and genetic factors related to EAT development and immunocyte infiltration. ceRNAs, competing for a shared pool of miRNAs, provide a complex network of interactions targeting the microRNA response elements (MREs) in protein-coding mRNAs, microRNAs, long noncoding RNAs, pseudogenic RNAs, and circular RNAs [25,26]. Recently, numerous studies have reported widespread functions, including RNA shearing, RNA editing, RNA modification and RNA binding, in human development and disease [25,26]. Thus, understanding the ceRNA networks linked to EAT development and immunocyte infiltration will lead to important insights into cardiac development c e d a b Fig. 4 The ceRNA network construction and detection related to infiltrating immunocytes. a The network plot presenting a regulatory landscape based on co-DE mRNAs, co-DE lncRNAs, and targeting miRNAs. b Fig. 3b showing Pearson's correlation coefficients for the hub mRNAs and candidate immunocytes. c, d After determination of the hub regulatory network, Fig. 3c shows a clustering analysis for hub mRNAs and lncRNAs, thus reflecting the sample discrimination ability. e Figure 3e demonstrates that the MIR210HG, CACNA2D4, and monocytes infiltration values were clustered  18:111 and immunoregulation involved in the age change of EAT in CHD patients in early life. Based on the complex bioinformatic analysis, our study demonstrated a framework to systematically characterize a regulatory map involving DE mRNAs, DE lncR-NAs, immunocyte infiltration, and potentially regulated pathways. In the comparisons among the neonate, infant, and child EAT samples from CHD patients, a total of 13 co-DE mRNAs and 47 co-DE lncRNAs were identified as core regulators that are primarily involved in the pathways of myoblast differentiation, ion transmembrane transport, and heart development. We subsequently constructed a ceRNA regulatory network for these co-DE mRNAs and co-DE lncRNAs as presented above. As a result, MIR210HG was identified as an important c d b a Fig. 5 The GSEA and GSVA in response to MIR210HG, CACNA2D4 gene expression, and monocyte infiltration. a, b After grouping, Fig. 4a shows the significant terms correlated with higher monocyte infiltration, higher CACNA2D4 expression, and higher MIR210HG expression, while Fig. 4b presents the opposite case. c, d Fig. 4c, d demonstrate the enriched genes of the PPAR signaling pathway and RIG-I-like receptor signaling pathway Page 10 of 13 Ma et al. J Transl Med (2020) 18:111 competing endogenous RNA of hsa-miR-449c-5p to regulate the CACNA2D4 expression level and thus regulate monocyte infiltration during the EAT development of CHD patients in early life. The PPAR signaling pathway was positively correlated with MIR210HG, CACNA2D4 expression, and monocyte infiltration, while the RIG-I like receptor signaling pathway was inversely correlated. Functionally, CACNA2D4 encodes a protein of the voltage-dependent calcium channel complex that serves as a mediator leading to the influx of calcium ions through the cytomembrane. Based on the PathCards database (https ://pathc ards.genec ards.org/) [27], we found that CACNA2D4 may be correlated with the pathways of T cell receptor signaling and dilated cardiomyopathy (DCM). Weeke et al. showed that the calcium channel subunit genes CACNB2 and CACNA2D4 were significantly correlated with familial atrial fibrillation [28]. While the regulation of cardiac rhythm has been studied, the exact mechanism has not been fully elucidated [28]. The long noncoding RNA-MIR210HG, as reported by recent studies, may function as an oncogenic regulator that promotes tumor metastasis, proliferation, and invasion in breast cancer, glioma, hepatocellular carcinoma, and non-small cell lung cancer [29][30][31][32][33][34]. In addition, Voellenkle's group and Lin's group demonstrated that MIR210HG may be involved in the pathway of hypoxic and inflammatory stress of endothelial cells [35,36]. Hsa-miR-449c-5p serves as an important mediator involved in brain development, gastric carcinoma growth, and liver cancer cell migration via extracellular signal and intracellular pathways that target the post-transcriptional regulation of genes [37][38][39]. Xu et al. found that the overexpression of hsa-miR-449c-5p may inhibit the osteogenic differentiation of valve interstitial cells via targeting Smad4 expression in calcific aortic valve disease [40].
The PPARs, including PPARα, PPARβ/δ, and PPARγ, were reported to regulate adipogenesis, lipid metabolism, and immunoinflammatory responses, especially in the maintenance of metabolic homeostasis [41]. Interestingly, recent findings have shown that PPARs, particularly PPARγ, are highly expressed in white adipose tissue (WAT) and BAT and mediate adipogenesis by regulated lipid oxidation, triglyceride accumulation, and glucose utilization under physiological and pathological conditions [41]. As the master regulator in adipocyte differentiation, PPARs also participate in angiogenesis, immune microenvironment construction  18:111 and energy metabolism in adipose tissue [41][42][43]. However, these processes of energy metabolism and construction of the immune microenvironment could be induced by EAT development and metabolic dysregulation due to disturbances in oxidative stress, ion transport, mitochondrial oxidative phosphorylation, and myocardial energy metabolism [44,45]. Wang et al. demonstrated that RIG-I may be a direct target and contribute to regulating the balance between the anti-inflammatory and proinflammatory responses and glucose homeostasis during the process of insulin resistance (IR), which is consistent with immune tolerance regulation [46]. In addition, the activation of the RIG-I like receptor signaling pathway may potentially link inflammation and cardiomyocyte apoptosis [47]. Renovato-Martins et al. demonstrate that monocyte infiltration into adipose tissue is a signature of organ early development, while macrophage infiltration is a hallmark of the inflammatory response in the adipose tissues of obesity patients [48]. Similarly, Wolf 's group found that the feedback loop between monocytes and macrophages not only provides immunological defense but also contributes to maintaining the function of EAT, which involves thermogenesis, metabolic balance, and local sympathetic nerve distribution [49]. Pirzgalska et al. demonstrated that monocytes-macrophages are intimately linked to the browning of white adipose tissue (WAT) and thermogenesis mediated by norepinephrine (NE) [50].

Conclusion
Related studies have proved that EAT development was regulated by the immune-inflammatory microenvironment and that the genetic regulation was complicated. In our study, these results further revealed the potential effects and research value of myoblast differentiation and ion transmembrane transport in EAT development during early childhood, and thus, these related DE genes were further assessed in the followup study. Additionally, a network including monocytes, MIR210HG, hsa-miR-449c-5p, and CACNA2D4 was detected in EAT development and may serve as a primary regulator controlling thermogenesis and immunoinflammatory response. These results might be used in novel strategies to develop effective therapies targeting myocardial energy supply and development in CHD patients. Our findings could further elucidate the molecular mechanisms involved in EAT, thereby enabling the exploration of EAT development and CHD pathophysiology from a new perspective, which may help improve the current therapies.