A Distinctive NAFLD Signature in Adipose Tissue from Women with Severe Obesity

Development and severity of nonalcoholic fatty liver disease (NAFLD) have been linked to obesity and white adipose tissue (WAT) dysfunction plays a key role in this relation. We compared the main features of subcutaneous (SAT) and visceral WAT (VAT) tissue dysfunction in 48 obese women without (Ob) and with NAFLD (Ob-NAFLD) undergoing bariatric surgery and matched for age, BMI and T2D status. Fat cell area, adipocyte size distribution, the degree of histological fibrosis and the mRNA expression of adipokines and genes implicated in inflammation, adipogenesis, angiogenesis, metabolism and extracellular matrix remodeling were measured by RT-qPCR in both fat depots. Ob-NAFLD group showed higher TG and lower HDL circulating levels, increased VAT fat cell area and similar WAT fibrosis in comparison with Ob group. A sPLS-DA was performed in order to identify the set of genes that better characterize the presence of NAFLD. Finally, we build a multinomial logistic model including seven genes that explained 100% of the variance in NAFLD and correctly predicted 100% of cases. Our data support the existence of distinctive NAFLD signatures in WAT from women with severe obesity. A better understanding of these pathways may help in future strategies for the prevention and treatment of NAFLD.


Introduction
NAFLD is defined as the presence of >5% of hepatic steatosis in the absence of other competing liver disease etiologies, concomitant use of alcohol or therapies that may induce liver steatosis [1]. NAFLD is the most prevalent chronic liver disease worldwide, affecting between 25-30% of the general population [2]. A "multiple hit" hypothesis has been proposed at the origin of its pathophysiology. Accordingly, multiple insults act together on genetically predisposed subjects to induce NAFLD [3]. These multiple parallel "hits" induce chronic systemic inflammatory processes altering key metabolic tissues, although complete portray of the factors that contribute to disease progression remain partially understood [4]. NAFLD patients tend to be obese and have other metabolic syndrome comorbidities. Thus, obesity is established as a contributor to the initial process leading to simple steatosis [5], indeed, more than 90% of severely obese patients that undergo bariatric surgery present with NAFLD [1]. However, NAFLD can also develop in subjects with normal body mass index (BMI) suggesting that metabolic alterations observed in obese patients that are found in other conditions might be drivers of NAFLD [6]. In fact, multiple studies have found that adipose tissue dysfunctions such as increased flux of free fatty acids (FFAs) to the liver, de novo hepatic lipogenesis and the release of pro-inflammatory cytokines by white adipose tissue (WAT) macrophages play principal roles in NAFLD development [7][8][9].
WAT expansion associated with body fat accrual leading to obesity, is accompanied by pro-inflammatory macrophage recruitment [10,11] which shifts cytokines production towards a more steatogenic, and fibrogenic profile. Indeed, it has been shown that the obesity-associated dysfunction in both subcutaneous (SAT) and visceral WAT (VAT) are responsible of up to 50% of circulating IL-6 secretion, contributing to systemic chronic inflammation [5]. Data on the importance of VAT in the chronic inflammation that drives NAFLD demonstrates that the lack of adipocyte AMP-activated protein kinase worsens liver disease by affecting brown and beige adipose tissue function [12].
WAT dysfunction is present in almost all patients with severe obesity and highly associated with NAFLD prevalence [13]. Nevertheless, few studies have systematically compared obesity-associated WAT dysfunction in relation to NAFLD condition. While early studies pointed to VAT depot as a main independent predictor of hepatic steatosis [14], liver inflammation and fibrosis [15], later investigations emphasized the association between NAFLD and SAT transcriptome [16][17][18] or macrophage phenotype [18]. Thus, du Plessis et al. built a predictive model to determine patients' liver histology based on SAT and VAT gene expression and found a model of 5 SAT genes as the most predictive in accurately assigning patients to histological groups. Results from this study underlined the importance of the SAT gene set at both early stages of NAFLD and in its progression to steatohepatitis [16]. Nevertheless, study groups were not matched for sex or central obesity. Recently, Song et al. constructed a transcriptomic signature in WAT which showed good performance in diagnosing NAFLD [19]. However, this model was based only on immunerelated genes and the depot origin of the WAT samples was not ascertained. Furthermore, none of these studies evaluated the degree of WAT fibrosis or fat cell size distribution. Therefore, we aimed to assess VAT and SAT histologic and transcriptomic differences of female patients with severe obesity and NAFLD compared to matched females for age, BMI and T2D status.

Results
The clinical characteristics of the 48 female participants are shown in Table 1. One-toone propensity score matching (PSM) methodology was used to obtain 24 obese women without (Ob) and 24 with NAFLD (Ob-NAFLD) pair-matched for age, BMI and T2D status, amongst subjects undergoing bariatric surgery between 2019 and 2021 at our Institution. After matching, serum triglycerides (TG) were higher and high-density lipoprotein (HDL) cholesterol levels were lower in the Ob-NAFLD group. In addition, two biochemical indices of liver fibrosis were calculated. Both FIB-4 and APRI (AST to platelet ratio index) scores were similar between groups and had low values, indicating great negative predictive value for advanced fibrosis.

Adipose Tissue Morphology
Histological analysis of adipocyte area and tissue fibrosis was performed on formalinfixed sections from both fat depots. Mean SAT-adipocyte area did not significantly differ among groups (Table 2, Figure 1A). On the contrary, mean VAT-adipocyte area was increased by 20% in the Ob-NAFLD group (Table 2, Figure 1A). Analysis of the adipocyte size distribution revealed that NAFLD condition is associated with a lower abundance (16% less) of smaller (<3000 µm 2 ) and higher abundance (55% more) of larger adipocytes (> 6000 µm 2 ) in VAT relative to Ob ( Figure 1B), while no significant differences were found in SAT. Furthermore, the mean ratio between individual SAT and VAT fat cell area was sig-nificantly lower in the Ob-NAFLD group (Table 2). Sirius red staining showed that fibrosis around adipocytes (i.e., pericellular fibrosis) were similar in both depots irrespective of NAFLD condition (Table 2, Supplementary Figure S1).

Differential Gene Expression Analysis
Descriptive statistics of the expression analysis of 75 genes implicated in WAT dysfunction evaluated both in SAT and VAT, grouped by NAFLD status are shown in Supplementary Table S1. The expression of 15 genes in SAT and eight genes in VAT was significantly different in Ob-NAFLD versus Ob subjects. Relative expression levels and fold changes between groups are reported in Table 3.

Differential Gene Expression Analysis
Descriptive statistics of the expression analysis of 75 genes implicated in WAT dysfunction evaluated both in SAT and VAT, grouped by NAFLD status are shown in Supplementary Table S1. The expression of 15 genes in SAT and eight genes in VAT was significantly different in Ob-NAFLD versus Ob subjects. Relative expression levels and fold changes between groups are reported in Table 3. In SAT, the pan-macrophage marker CD68 and the M2-type macrophage markers MSR1/CD204 and MRC1/CD206 were found upregulated in patients with NAFLD (Table 3, Supplementary Figure S2). The expression of vascular endothelial growth factor B (VEGFB) and angiopoietin 1 (ANGPT1), two pro-angiogenic factors was also increased in the Ob-NAFLD group by 38% and 41%, respectively. In addition, leptin (LEP) gene expression was up-while its receptor (LEPR) was downregulated. Several genes implicated in ECM composition and remodeling were found altered in the NAFLD group. Thus, despite similar levels of histological fibrosis, TGFB1 gene expression was upregulated by 70%. Similarly, the expression of the α-1 chains of collagens I, IV and VI were up-modulated, especially COL6A1 (by 940%). On the contrary, the tissue inhibitor of metalloproteinases-3 (TIMP3) was down-regulated by 48%.
In VAT, the expression of hypoxia-inducible factor 1-alpha (HIF1A), ANGPT2 and the thermoregulatory gene cell death-inducing DFFA-like effector A (CIDEA) was found slightly decreased in NAFLD obese females (Table 3, Supplementary Figure S3). The expression of senescence genes TP53/P53 and CDKN2A/P16 was decreased (by 49 and 35%, respectively) while LEP and monoglyceride lipase (MGLL) were significantly upregulated (by 87 and 67%, respectively). Table 3. Significant results of differential gene expression analysis in SAT and VAT from Ob-NAFLD vs. Ob groups.

Principal Component Analysis (PCA) to Identify NAFLD Patients
The PCA including expression data of all 150 genes did not spot defined cluster differences between patients with and without NAFLD (Q2X = 21.9%5, R2X = 45.3%) (Supplementary Figure S4). In order to find an adipose tissue-associated genomic signature capable to differentiate patients with NAFLD, a sPLS-DA model cross-validated via LOO (R2Y = 0.94; Q2 = 0.43, p = < 0.001) (Figure 2A Table S1). All genes fro the signature model had mean expressions significantly different in patients with NAFL Multinomial logistic regressions were performed to assess the effects of SAT an VAT gene expression on the likelihood that patients have NAFLD (Table 4). SAT-TGF and VAT-P53 resulted the only genes independently associated with the presence NAFLD. Patients with higher expressions of SAT-TGFB1 were 4.04 (1.15-14-23) tim Multinomial logistic regressions were performed to assess the effects of SAT and VAT gene expression on the likelihood that patients have NAFLD (Table 4). SAT-TGFB1 and VAT-P53 resulted the only genes independently associated with the presence of NAFLD. Patients with higher expressions of SAT-TGFB1 were 4.04 (1.15-14-23) times more likely to present NAFLD whereas increasing VAT-P53 expression was associated with a 0.25 likelihood reduction (0.07-0.85) of having the disease. A model including these genes explained 45.2% of the variance in NAFLD and correctly classified 80.0% of cases (p = 0.000). Moreover, we found a second model including 7 genes ( Table 4) that explained 100% of the NAFLD variance (p = 0.000) though no gene reached independent significance explained because of gene intercorrelation (Supplementary Table S2).

Discussion
In the present study, we characterized the main histological and gene expression differences in SAT and VAT between females with severe obesity without and with NAFLD as assessed from hepatic ultrasound. Through matching patients by age, BMI and T2D, WAT alterations induced by NAFLD per se could be evaluated avoiding potential confounders. Overall, anatomical changes were mainly found in VAT adipocyte area of NAFLD patients with a shift in the proportion of adipocytes to a larger fat cell area range. However, in transcriptomic analysis, both fat compartments had differentially expressed genes associated with NAFLD, which constituted a highly discriminant WAT-signature of 25 genes. Finally, SAT-TGFB1 and VAT-P53 expression were genes with high likelihood estimates independently associated with NAFLD.
Clinically, obesity associates increased morbidity and mortality when combined with NAFLD. Pathogenetically, obesity and the ensuing insulin resistance contribute to the initial fat accumulation in the hepatocyte and to its progression into nonalcoholic steatohepatitis (NASH), NASH-related cirrhosis and hepatocellular carcinoma [20]. Nevertheless, obesity does not always concur with NAFLD, as not all NAFLD patients are obese and not all subjects living with obesity present NAFLD [21,22]. Different causal links have been described to explain the association between obesity and NAFLD [23][24][25], and altered adipose tissue biology has been recognized as a key early event in its initiation. However, to our knowledge, the number of studies systematically comparing the biology of WAT in relation to the presence of NAFLD is very limited to date [16,26,27], comparisons are made with respect to healthy lean controls [28] or are focused on specific WAT-derived circulating factors [29,30].
Obesity-induced changes in adipokine secretion can modify the adipose tissue-liver crosstalk [31,32] and the ratio between leptin and adiponectin has been proposed as a superior biomarker of NAFLD than adiponectin or leptin measurements alone [33,34]. Here we report higher SAT and VAT leptin expression in NAFLD group while comparable levels of adiponectin expression, what translates into increased visceral leptin-to-adiponectin ratio in this group. Greater leptin expression in both fat depots contributed to classify patients with NAFLD in the sPLS-DA, although it was not included in the multinomial logistic regression model.
The 'portal theory' suggested that, in obese patients, venous drainage of increasing amounts of pro-inflammatory factors and free fatty acids from VAT to the liver via the portal system favors the development of hepatic insulin resistance and liver steatosis [35][36][37].
However, some authors presented evidence challenging this view and questioned its relative contribution with respect to SAT [38]. Despite VAT dysfunction is generally considered as the main contributor to NAFLD [14,15], SAT expression of genes was overrepresented in all group comparisons in our study, being 15 vs. 8 the number of genes in SAT vs. VAT, respectively, whose expression was significantly different in the two sample comparisons; 18 vs. 7 were included in the sPLS-DA model and 4 vs. 3 were included in the multinomial logistic regression model. However, this does not imply a direct causal relationship of SAT dysfunction in the development of NAFLD and the design of this study prevents drawing mechanistic conclusions. One possibility is that greater SAT dysfunction leads to greater VAT dysfunction, and this in turn is associated with the development of the NAFLD.
In a multicenter study, du Plessis et al. built a predictive model of NAFLD based on SAT and VAT gene expression levels from severely obese patients [16]. Interestingly, their results emphasize the significance of SAT inflammation and the highly predictive nature of the SAT gene set in classifying liver histology over VAT. In this sense, previous studies described that SAT macrophages resembled the pro-inflammatory phenotype of VAT macrophages and were significantly increased in patients with NASH and fibrosis [18]. Similarly, the transcription of pro-inflammatory genes and macrophage numbers in SAT correlated with hepatic fat content in other study [17].
Our model including SAT expression of TGFB1 and VAT expression of P53, two tumor suppressor genes, correctly classified 80.0% of NAFLD cases. TGFB1, which was upregulated in Ob-NAFLD group, promotes the release of inflammation mediators, remodeling and collagen deposition in WAT [39][40][41]. TGFB1 release from SAT is increased in obesity [42] while systemic blockade of its signaling protects mice from obesity, diabetes and hepatic steatosis [43].
Accumulating evidence suggests that the role of P53 in the progression of NAFLD is complex and context-dependent [44]. P53 plays multiple functions in cell cycle, cellular senescence, autophagy, apoptosis, and metabolism [45]. While hepatic p53 expression is elevated in patients with NASH [46,47] and positively correlates with the degree of steatosis [47], WAT p53 seems to play a specific role as a negative regulator of adipogenesis [45] and reduces lipid droplet accumulation [48]. In our study, visceral expression of P53 was reduced in NAFLD patients but, strikingly, this was accompanied by a higher adipocyte hypertrophy.
Interestingly, the expression of CIDEA, which was found downregulated in VAT and main predictor of NAFLD, has been positively associated with human obesity but inversely related to NAFLD severity [49]. CIDEA is part of the cell death-inducing like effector (CIDE) family that regulates hepatic lipid homeostasis controlling lipid droplet growth. The downregulation of VAT-CIDEA found in our NAFLD population points out one of the genes that might be independently associated to liver injury irrespective from obesity.
SAT-UCP2 which was highly increased in patients with NAFLD and was the third gene in the signature contributors, is a mitochondrial anion carrier which expression has been associated to cells that are exposed to obesity-associated oxidative stress as a defensive mechanism [50]. UCP2 overexpression has been found to induce acute liver injury and HFD-induced liver damage while its downregulation its associated with NAFLD amelioration in animal models [51].
MGLL, whose visceral expression was increased in patients with NAFLD, its part of the monoglyceride lipases family that hydrolyzes triglycerides to fatty acids and glycerol in VAT [52] and may indicate an increased lipolysis in the NAFLD group, leading to excessive circulating FFA levels and contributing to the development of NAFLD.
We report a decreased expression of the autophagy related gene 5 (ATG5) in SAT from NAFLD patients. While pharmacological inhibition, silencing or knockdown of ATG5 in hepatocytes results in increased triglyceride levels and lipid droplet accumulation [53], WAT expression of ATG5 is increased in obesity, especially in the VAT depot [54], and is a key regulator of adipogenesis. In mature adipocytes, however, is still not clear the role of autophagy and if it is beneficial or detrimental [55].
In contrast with our results of VAT-ANGPT2 being downregulated in patients with NAFLD, a study of 93 severely obese subjects found that ANGPT2 expression in VAT was associated with different NAFLD features, including steatosis, ballooning, portal and lobular inflammation [56]. However, this correlation was lost after correction for T2D and insulin resistance and plasma ANGPT2 levels were not related to NAFLD features. Moreover, the study did not include a control group without NAFLD.
Finally, we found an upregulation of mannose receptor C-type 1 (MRC1) in SAT from the NAFLD group. MRC1, a surface marker of M2 or anti-inflammatory macrophages, was previously found increased in SAT from obese subjects with NAFLD versus obese with normal intrahepatic triglyceride content and lean control patients [57].
We acknowledge our study is not without limitations. First, diagnosis of NAFLD was based only on ultrasound assessment, and not in the gold standard liver biopsy. Thus, presence of mild liver steatosis could not be ruled out in our study subjects classified in the non-NAFLD group. Although liver biopsies were not performed, data from validated scores would suggest our study patients were at early stages of NAFLD. Second, even though the groups were paired for BMI and T2D, body composition and fat distribution of the participants were not evaluated. Finally, the approach presented here, despite highly discriminant to diagnose NAFLD, would not be a more feasible alternative than the gold standard liver biopsy nor than non-invasive tests such as biomarkers or ultrasound, since a VAT sample would be still needed.
Overall, within its limitations, our data support the existence of a distinctive NAFLDassociated histological and transcriptomic signature in WAT from females with severe obesity even when compared to patients with similar clinical parameters. Further studies are needed to confirm our data, as well as functional studies are warranted to establish causal relationships between these genes and the development of NAFLD. A better understanding of this pathways may help in future strategies for the prevention and treatment of NAFLD.

Study Population
Forty-eight severely obese women were selected among candidates undergoing bariatric surgery at the Obesity Unit of the Hospital Clínic of Barcelona. The study's exclusion criteria were a history of malignancy, chronic inflammatory diseases, active infectious diseases, drug abuse or daily alcohol consumption >20 g. Twelve patients underwent a laparoscopic sleeve gastrectomy and thirty-six underwent a gastric bypass. Two matched female cohorts, matching NAFLD 1:1 to non NALFD subjects (24 vs. 24), were obtained using "nearest neighbor" matching for age, BMI and T2D and the maximum allowed distance was a ∆ of 0.001 (SPSS, version 25.0, Chicago, IL, USA). The matching significantly reduced differences in these variables. Ethics committee approval conforming to the Declaration of Helsinki was obtained from the Clinical Research Ethics Committee (CEIC) of Hospital Clinic de Barcelona. All participants provided written informed consent.

Clinical and Anthropometric Data
The patients' anthropometric measurements were collected in the same consultations following standardized procedures and hematological and biochemical parameters were determined at the Core Laboratory of the Biomedical Diagnostic Center using an Advia 2400 analyzer (Siemens Healthcare S.L.U., Getafe, Spain). Anthropometric and clinical data are summarized in Table 1.

Determination of NAFLD
Presence of NAFLD was ascertained as part of the pre-bariatric surgery evaluation and following exclusion of secondary causes of liver steatosis and excessive alcohol intake (defined as ≥20 g/d) [58]. NAFLD was defined as the presence of hepatic steatosis on ultrasonography (US). Hepatic US was performed in all patients after 6 h fasting, by a single experienced using a clinical US scanner (Aplio i-800, Canon Medical Systems S.A., Madrid, Spain) equipped with a i8CX1 1-8-MHz curved US transducer used for conventional B-mode examination. Each subject was examined in the supine and left lateral positions during quiet inspiration. Hepatic steatosis was diagnosed on the basis of characteristic US features: evidence of diffuse hyperechogenicity of the liver relative to the kidneys, ultrasound beam attenuation and poor visualization of intra-hepatic vessel borders and diaphragm. Semiquantitative US scoring for the degree of hepatic steatosis was not available in this study. Absence of advanced liver fibrosis was evaluated using the fibrosis-4 (FIB-4) and the AST to platelet ratio index (APRI) according to the proposed cut-offs [58].

White Adipose Tissue Biopsies
Paired SAT and VAT samples were obtained at the time of surgery. Samples were collected in DMEM and rinsed in PBS. A portion was immediately frozen before RNA analysis. Other portion was fixed overnight at 4 • C in 4% paraformaldehyde and processed for standard paraffin embedding. Starting at the tissue apex 3 × 3 µm sections were made at a minimum of 100 µm intervals across the sample tissue. Serial sections were matched for additional independent analyses.

Morphometry and Histopathology
Hematoxylin and eosin (H&E) staining was used to assess adipocyte morphology. Digital images were captured under an X600 microscope (Olympus, Tokyo, Japan) at 4× magnification. Adipocyte size of at least 3000 cells per sample was measured within a minimum of 5 micrographs from randomly selected fields using Adipocytes Tools, an ImageJ macro-based algorithm for ImageJ software (National Institutes of Health, Bethesda, MD, USA; http://imagej.nih.gov/ij/, last accessed 20 August 2021). Adipocyte average area was calculated and frequency distribution analysis into bin intervals of 200 µm 2 was performed.
Sirius red staining was used for quantification of pericellular fibrosis (i.e., extracellular matrix accumulation around the cells). Automated analysis of at least 10 images per sample at 10× magnification has been carried out using MRI Fibrosis Tool, an ImageJ macro-based algorithm, and expressed as a percentage of red staining (fibrosis)/tissue surface ratio. SAT:VAT fibrosis ratio has been calculated as the ratio between the percentage of fibrotic area in SAT and VAT from each patient.

RNA Extraction and Real Time PCR
Total RNA was isolated using RNeasy Lipid Tissue Mini Kit (Qiagen, Hilden, Germany). Concentration and purity were measured using a NanoDrop 1000 spectrophotometer (Thermo Scientific, Waltham, MA, USA). Equal amounts of RNA from SAT and VAT (2 µg) were reverse-transcribed using the Superscript III RT kit and random hexamer primers (Invitrogen, Carlsbad, CA, USA). Reverse transcription reaction was carried out for 90 min at 50 • C and an additional 10 min at 55 • C. An expression analysis of 75 genes involved in WAT dysfunction and related to inflammation, adipogenesis, autophagy, fatty acid metabolism and oxidation, adipocyte brightening, glucose metabolism and adipokines was performed in both fat depots. Real time quantitative PCR (qPCR) was performed with a 7900HT Fast Real-Time PCR System (Applied Biosystems, Foster City, CA, USA) using GoTaq ® qPCR Master Mix (Promega Biotech Ibérica, Madrid, Spain). Expression relative to the housekeeping gene RPL6 was calculated using the delta Ct (DCt) method. Gene expression is presented as the 2ˆ(-DCt) values. Fold changes were calculated as the ratio between the mean expression levels in Ob-NAFLD subjects and Ob subjects. The list of primers used in this study is provided in Supplementary Table S4.

Statistics
Continuous data with normal and non-normal distribution is expressed with arithmetic means and standard deviations (SD), or medians and interquartile ranges (IQR), respectively. Categorical variables are expressed with frequencies and proportions. No data transformation was performed on gene expression when performing univariate analysis and normality assumption was tested with Shapiro-Wilk test. Linearity, and absence of multicollinearity were also checked. For each gene, the magnitude of the expressiondifference between groups was calculated with mean fold changes (FC) and tested with Mann-Whitney U test, Welch's t-test or Student's t-test when adequate.
Multivariate dimensionality reduction methods were conducted to build a predictive classification model for NAFLD. First, missing gene expression values (n = 886, 12.02%) were imputed with the k-nearest neighbor algorithm and standardized to zero means and unit variances giving similar importance to all genes independently to its intrinsic tissue expression. Principal component analysis (PCA) using singular value decomposition was implemented to examine the intrinsic dimension and visualize the general structure of the transcriptome across groups. To find the most suitable signature, leave-one-out cross-validation of a sparse partial least square discriminant analysis (PLS-DA) was used to determine the optimal number of components and genes to be included in the model. The number of components was chosen based on the estimation of the lowest balanced error rate (BER) and the genes based on the lowest prediction error for each subset of genes on each component. R2Y (the sum of squares) and Q2Y (the predictive performance) values were assessed to ensure the absence of overfitting of the final model. Area under the curve (AUC) values were calculated from the predicted scores in the LOO cross-validation process minimizing the risk of overfitting. A clustered image map was created employing multivariate Euclidean distance metric with complete linkage method and presented with associated dendrograms. Finally, Spearman correlations followed by polynomial regressions with enter method were then performed to identify genes independently associated with the presence of NAFLD among those included in the model. Correction by age adjustment was used in all multivariable models.
All the comparisons stated as different in the present manuscript have statistical significance with a two-tailed p-value < 0.05. GraphPad PRISM (GraphPad Software, version 6.0, San Diego, CA, USA), Statistical Package for Social Sciences software (SPSS, version 25.0, Chicago, IL, USA) and R (The R Project for Statistical Computing, version 4.1, Vienna, Austria) software environment [59] were used to perform the analyses.  Informed Consent Statement: Informed consent was obtained from all subjects involved in the study.

Data Availability Statement:
All data presented in this study are reported in this manuscript or available in Supplementary Material.