Endothelial-immune crosstalk contributes to vasculopathy in nonalcoholic fatty liver disease.

The top cause of mortality in patients with nonalcoholic fatty liver disease (NAFLD) is cardiovascular complications. However, mechanisms of NAFLD-associated vasculopathy remain understudied. Here, we show that blood outgrowth endothelial cells (BOECs) from NAFLD subjects exhibit global transcriptional upregulation of chemokines and human leukocyte antigens. In mouse models of diet-induced NAFLD, we confirm heightened endothelial expressions of CXCL12 in the aortas and the liver vasculatures, and increased retention of infiltrated leukocytes within the vessel walls. To elucidate endothelial-immune crosstalk, we performed immunoprofiling by single-cell analysis, uncovering T cell intensification in NAFLD patients. Functionally, treatment with a CXCL12-neutralizing antibody is effective at moderating the enhanced chemotactic effect of NAFLD BOECs in recruiting CD8+ T lymphocytes. Interference with the CXCL12-CXCR4 axis using a CXCR4 antagonist also averts the impact of immune cell transendothelial migration and restores endothelial barrier integrity. Clinically, we detect threefold more circulating damaged endothelial cells in NAFLD patients than in healthy controls. Our work provides insight into the modulation of interactions with effector immune cells to mitigate endothelial injury in NAFLD.

Introduction

Nonalcoholic fatty liver disease is the most common liver disease in developed countries, affecting up to a third of western populations, and up to 40% in South East Asia (Muthiah & Sanyal, 2020). The spectrum of NAFLD includes the more benign nonalcoholic fatty liver, also known as simple steatosis, and the more progressive form of nonalcoholic steatohepatitis (NASH) and can progress to fibrosis (Chalasani et al, 2018). Nonalcoholic steatohepatitis is characterized by lobular inflammation and ballooning and can progress to develop fibrosis. While NAFLD remains a liver disease, the top cause of mortality in patients with NAFLD is cardiovascular mortality, which occurs independently of the shared metabolic risk factors such as insulin resistance and obesity (Targher et al, 2020). Patients with NAFLD are at increased risk of developing multiple vascular complications, including coronary artery disease, cerebrovascular disease, and peripheral vascular disease.

There is emerging evidence of the link between NAFLD and impaired vascular health. Clinical measures of vascular function have established the association of reduced flow‐mediated vasodilatation, altered carotid artery intimal medial thickness, coronary calcification, and low coronary flow reserve with NAFLD (Sookoian & Pirola, 2008; Federico et al, 2016). In line with this, regression of NAFLD is associated with a lower risk of carotid atherosclerosis in observational cohort studies (Sinn et al, 2016, 2017). Even in the absence of major cardiometabolic risk factors, healthy individuals demonstrated better vascular functions than NAFLD patients (Long et al, 2015; Al‐Hamoudi et al, 2020), supporting that NAFLD alone could contribute to the dysfunction of systemic vasculatures. Various mechanisms have been proposed as plausible evidence for accelerated atherosclerosis and increased cardiovascular risks in NAFLD patients, including a high oxidative stress state, macrophage activation (Targher et al, 2010), and systemic release of proatherogenic and thrombotic factors like tumor necrosis factor‐α, interleukin‐6, and oxidized low‐density lipoprotein (Stols‐Goncalves et al, 2019). Many of these paracrine mediators are known to activate endothelial cells and potentially result in endothelial dysfunction. Inflammation may also augment lipid risks to further drive atherosclerosis (Libby, 2021). However, there remain knowledge gaps in NAFLD‐related endothelial disease biology that are not resolved by the current evaluation of vascular function using imaging modalities and serum biomarkers.

We aim to understand the molecular basis of endothelial pathophysiology in NAFLD‐associated vasculopathy. Due to the difficulties of obtaining fresh vascular tissues from patients, we harnessed the use of blood outgrowth endothelial cells (BOECs) to develop personalized endothelial models from NAFLD patients and healthy subjects. BOECs originate from endothelial colony‐forming cells, which are bone marrow–derived progenitors found in the circulation and within vascular endothelium (Yoon et al, 2005). They are distinct from the early endothelial progenitor cells, which form a heterogeneous culture of monocyte‐derived cells from hematopoietic stem cell origin, whereas BOECs give rise to a homogeneous population of mature endothelial cells with cobblestone morphology and are capable of tube formation. While the precursor of BOECs, endothelial colony‐forming cells, have been explored for therapeutic revascularization, BOECs prove to be a useful surrogate of patient endothelial cells to investigate the biology of vasculopathy in various disease contexts such as pulmonary arterial hypertension, diabetes, and ischemic heart disease (Paschalaki & Randi, 2018). Patient‐derived BOECs were able to recapitulate clinical phenotypic alterations, capture the complexities of genetics and/or epigenetics, and demonstrate differential responses to disease‐relevant stressors (Hebbel, 2017). Here, we utilized BOECs to facilitate experimentations including transcriptomic analysis, mechanistic interrogation, and endothelial‐immune cell coculture assays. We also analyzed circulating endothelial cells (CECs), which were shed from endothelial lining into the bloodstream following vascular damage (Hebbel, 2017), to yield insights on endothelial injury in NAFLD.

Results

Development of human blood outgrowth endothelial cells for disease modeling

Patients with NAFLD were diagnosed with at least 5% steatosis characterized on liver biopsies, while selected healthy controls had a Controlled Attenuation Parameter score of less than 248 on vibration controlled transient elastography, indicative of less than 5% steatosis (more details in Materials and Methods under “patient enrollment”). We employed established protocols (Martin‐Ramirez et al, 2012b; Ormiston et al, 2015) to develop BOECs from NAFLD patients and healthy subjects. Briefly, peripheral blood mononuclear cells (PBMCs) isolated from whole blood samples were cultivated to derive BOECs (Fig 1A). We then characterized the growth dynamics, marker expressions, and functions of NAFLD and healthy control BOECs. BOEC colonies usually emerged after 7–14 days of PBMC cultivation (Fig 1B). The colonies were isolated, passaged into endothelial growth media, and further cultured to stabilize the BOEC lines that displayed cobblestone endothelial cell morphology. Both NAFLD and healthy BOECs had comparable doubling times (Fig 1C). More than 90% of BOECs expressed endothelial markers, including CD31 (PECAM1), CD144 (CDH5), and CD146, but negligible expressions for leukocyte markers CD45, CD14, and CD68 and progenitor cell marker CD133 (Fig 1D), suggesting high purity of our BOEC derivation. Immunostaining confirmed their positive expressions of endothelial markers (Fig 1E). Functionally, the BOECs were able to undergo angiogenesis, typical of endothelial property, in three‐dimensional fibrin gel bead‐based sprouting assays. Sprouts could be observed as early as 6 h (Fig 1F, representative images). By 24 h, both NAFLD and healthy BOECs showed a comparable number of sprouts per bead and sprout lengths (Fig 1F, lower panel). During sprouting angiogenesis, endothelial cells undergo a series of dynamic changes in their tip‐ and stalk‐like cell phenotypes to facilitate cell migration, proliferation, and stabilization, which ultimately leads to sprout formation (Blanco & Gerhardt, 2013). Phenotypically, we did not observe differences in the angiogenic capacities between NAFLD and healthy BOECs. However, at the molecular level, sprouted NAFLD BOECs might be primed for angiogenesis as demonstrated by their significantly higher gene expressions for tip cell markers than sprouted healthy BOECs (Fig EV1). Concomitantly, the role of pathological angiogenesis in NAFLD is well known, likely driven by pro‐angiogenic factors arising from chronic inflammation, tissue hypoxia, and hepatocyte injury (Coulon et al, 2012; Hammoutene & Rautou, 2019; Lefere et al, 2019). Taken together, these quality control characterizations in Figs 1 and EV1 confirmed mature endothelial phenotypes of both NAFLD and healthy BOECs. Good quality cell lines were used for downstream experimentations.

Figure 1

Generation and characterization of blood outgrowth endothelial cells from healthy and NAFLD subjects

Workflow schematic of BOEC generation from PBMC samples isolated from NAFLD patients and healthy donors.

Upper panel: Representative images of BOEC colonies from PBMC cultivation 10 days post‐seeding. Dotted line in each image outlines a BOEC colony. Lower panel: Stabilized BOEC cultures after passaging of colonies (scale bar, 250 µm).

Proliferation doubling time for NAFLD and healthy BOECs based on passages 2 and 3. Sample sizes are n = 9 healthy, n = 7 NAFLD; with 3 biological replicates per donor. Box‐whisker plots indicate median (middle line), 25th, 75th percentile (box), and the lowest/highest data points (whiskers). ns—not significant (t‐test).

Left panel: Flow cytometry characterization of BOECs for endothelial, leukocyte, and progenitor cell markers (gray—isotype control; red/blue/green/purple—cell lineage marker staining). HUVEC and monocytic THP‐1 cells were used as positive controls for endothelial and leukocytic markers, respectively. Right panel: Percentages of positively expressing cells. Sample sizes are n = 4–6 healthy, n = 2–7 NAFLD. Results are indicated by mean ± SD. *P < 0.05; **P < 0.01; ****P < 0.0001; ns—not significant (one‐way ANOVA).

Representative images of immunostaining for endothelial markers, PECAM1 and CDH5 (scale bar, 100 µm).

Left panel: Longitudinal monitoring of BOEC angiogenesis in fibrin gel bead‐based sprouting assays over 24 h (representative images; scale bar, 100 µm). Right panel: Box‐whisker plots of quantifications of sprout lengths and numbers of sprout per bead at 24 h. n = 3 healthy, n = 3 NAFLD, with 4–5 beads analyzed per donor and 5–11 sprouts per bead. ns—not significant (t‐test). Box‐whisker plots indicate median (middle line), 25th, 75th percentile (box), and the lowest/highest data points (whiskers).

Figure EV1

Gene expressions of tip‐stalk markers in sprouting BOECs

Expressions were determined by qPCR and normalized to monolayers (nonsprouting). Data presented as mean ± SD. n = 3 donors; ns—not significant; ****P < 0.0001 (t‐test).

NAFLD endothelial cells show enhanced chemokine hallmarks

To explore the molecular basis of endothelial pathophysiological mechanisms in NAFLD, we performed RNA‐sequencing on BOECs derived from 3 NAFLD patients (2 females and 1 male, 37.7 ± 11.5 years old) and 3 healthy controls (2 females and 1 male, 41.7 ± 8.6 years old) (Appendix Table S1). Principal component analysis illustrated a clear segregation of transcriptomic profiles between NAFLD and healthy groups (Fig 2A). A total of 670 differentially expressed genes were identified, including 535 upregulated genes and 135 downregulated genes in NAFLD BOECs in comparison to healthy BOECs (Fig 2B; false discovery rate < 0.05, P < 0.05; fold change > 2 or < −2, respectively). Gene ontology analysis revealed that NAFLD BOECs are largely enriched in processes related to cell locomotion, extracellular matrix organization and chemotaxis, whereas filament organization and mitogen‐activated protein kinase (MAPK) cascade were most prominent in healthy BOECs (Fig 2C). Based on the upregulated genes in NAFLD BOECs, chemotaxis was identified as the most central network (Cluster 1) by molecular complex detection (MCODE), interconnecting with cell chemotaxis (leukocyte migration), cytoskeleton organization, and cell morphogenesis (Fig 2D). In contrast, no MCODE cluster was found using the downregulated gene set (Fig EV2A), suggesting that the upregulated genes in NAFLD BOECs might exert a stronger biological effect due to intertwined processes. Correspondingly, heatmap visualization of the genes from the top enriched network showed that NAFLD BOECs expressed higher levels of chemokines including the CC, CXC, and CX3C families than healthy BOECs (Fig 2E).

Figure 2

Transcriptomic analysis of NAFLD and healthy BOECs

Principal component analysis of NAFLD (n = 3) and healthy (n = 3) BOEC transcriptomes, with 3 biological replicates per donor.

Volcano plot depicting differentially expressed genes comparing NAFLD BOECs with healthy BOECs. Red and blue dots represent upregulated and downregulated genes in NAFLD BOECs, respectively.

Enriched gene ontology terms of the differentially expressed genes.

Metascape enrichment network visualization showing intracluster and intercluster interactomes of the upregulated genes in NAFLD BOECs. Nodes that share the same cluster annotations are close to each other. Densely connected clusters were identified using MCODE algorithm and outlined in red.

Heatmap featuring the genes in MCODE Cluster 1.

STRING protein–protein interaction network were analyzed by maximum neighborhood component, showing hub genes as ranked and indicated by gradient of red colors. Blue nodes represent differentially expressed genes in the extended networks.

Figure EV2

Transcriptomes of healthy and NAFLD BOECs

Metascape enrichment network visualization of interactomes of the downregulated terms in NAFLD BOECs. No identifiable MCODE complexes were found.

Ten MCODE complexes identified in the upregulated networks in NAFLD BOECs, ranked, and colored by their identities. The identities of the complexes were listed in the table.

To prioritize candidate genes for further mechanistic interrogation, we used 12 topological data analysis methods to identify and rank top 20 genes per method, which might play a central role in the network interactome (Appendix Table S2). Chemokines including CXCL3, CXCL5, CXCL6, CXCL10, CXCL11, CXCL12, CCL20, and CX3CL1 were consistently reflected by at least 6 topological analysis methods. Of those, CXCL10, CXCL12, and CX3CL1 were ranked as top 3 genes by at least three methods, emphasizing their importance in the network. Fig 2F is a visualization of one of the methods used. We also applied the MCODE algorithm that verified the “chemokine‐mediated signaling pathway” as the top MCODE complex containing most of the aforementioned chemokine genes (Fig EV2B). Taken together, enhanced chemokine signatures could be a hallmark for NAFLD endothelial cells.

NAFLD patient endothelial cells and vascular endothelia of in vivo disease models demonstrate intensified CXCL12 levels

We selected several chemokines from the prioritized list of candidate genes for further validation in a larger number of endothelial cell lines. Gene expression profiling demonstrated that CXCL10 and CXCL12 were significantly upregulated in NAFLD BOECs (n = 10) compared with healthy BOECs (n = 12), while the other chemokines were not significantly different at baseline (Fig 3A, donor demographics in Appendix Table S3). To examine these chemokine expressions under a relevant pathophysiological context, we introduced autologous plasma as a stressor to the BOEC cultures. Upon exposure of BOECs to their autologous plasma, similarly, CXC chemokine ligands were significantly upregulated in NAFLD BOECs than in healthy BOECs (Fig EV3A). It was reported that levels of free fatty acids (FFAs) and lipopolysaccharide (LPS) in plasma were positively correlated with NAFLD (Zhang et al, 2014; Zhu et al, 2015). These factors are known to activate endothelial cells and exacerbate inflammation. We further examined the chemokine expression profile of BOECs in the presence of FFAs and LPS. There were no significant differences in lipid uptake by both NAFLD and healthy BOECs (Fig 3B). Under exposure to FFAs and LPS, NAFLD BOECs had visibly more pronounced CXCL12 and CX3CL1 expressions than healthy controls (Fig 3C). NAFLD BOECs seemed intrinsically primed with higher CXCL12 expression than healthy BOECs.

Figure 3

Validation of chemokine hallmarks in human endothelial cells and mouse disease models

Fold change in gene expressions of CXCL10, CXCL11, CXCL12, CCL20, and CX3CL1 in NAFLD (n = 10) BOECs, normalized to healthy (n = 12) BOECs, with 2 technical replicates per donor. Results are indicated by mean ± SD. *P < 0.05; **P < 0.01; ns—not significant (Mann–Whitney test).

Left: Representative images showing nontreated and LPS + FFA‐treated NAFLD and healthy BOECs, stained with Nile Red (scale bars, 100 µm). Right: Fluorescence intensity of intracellular Nile Red stain was quantified at 515/585 nm, and readings were normalized with blanks (n = 3 donors/group). Results are indicated by mean ± SD. ns—not significant (two‐way ANOVA comparing NAFLD and healthy groups).

Fold change in gene expressions of CXCL10, CXCL11, CXCL12, CCL20, and CX3CL1 in LPS + FFA‐treated or nontreated NAFLD (n = 5) and healthy (n = 5) BOECs, normalized to nontreated healthy BOECs, with 3 technical replicates per donor. Results are indicated by mean ± SD. *P < 0.05; **P < 0.01; ***P < 0.001; ****P < 0.0001 (two‐way ANOVA).

Protein concentrations of chemokines in the conditioned culture media of healthy and NAFLD BOECs. *P < 0.05; ns—not significant (Mann–Whitney test). Box‐whisker plots indicate median (middle line), 25th, 75th percentile (box), and the lowest/highest data points (whiskers).

Left: Representative images (Bright field; DAPI, nuclei; CDH5, endothelial cells; CXCL12) of liver sections from chow diet‐fed and HFHC diet‐fed humanized mice (scale bar, 40 μm). Boxed regions were further magnified (scale bar, 30 µm). Right: Box‐whisker plot of Mander’s coefficient (M2) that was a measure of the degree of CXCL12 overlap with CDH5‐positive vasculatures in the livers. Box‐whisker plots indicate median (middle line), 25th, 75th percentile (box), and the lowest/highest data points (whiskers). Chow diet‐fed (n = 4) and HFHC diet‐fed (n = 4) humanized mice were analyzed, with 10–20 independent regions of interest per mouse were analyzed for image quantification. ****P < 0.0001 (Mann–Whitney test).

Left: Representative images (Bright field; DAPI, nuclei; CDH5, endothelial cells; CXCL12) of aortic sections from chow diet‐fed and LIDPAD mice (scale bar, 40 μm). Boxed regions were further magnified (scale bar, 30 µm). Arrows point at aortic endothelial cells. AW—aortic wall; L—lumen. Right: Box‐whisker plot of CXCL12 fluorescent intensity in the aortic vascular endothelia. CXCL12 fluorescent intensity for each delimitated endothelial region is represented as individual data points on the plot. Box‐whisker plots indicate median (middle line), 25th, 75th percentile (box), and the lowest/highest data points (whiskers). Chow diet‐fed (n = 3) and LIDPAD diet‐fed (n = 3) mice were analyzed, with 10–20 delimitated endothelial regions per mouse analyzed for image quantification. ****P < 0.0001 (Mann–Whitney test).

Left: Representative images showing aortic sections harvested from mice fed with either chow or LIDPAD diet. Sections were stained for leukocyte markers CD45 (red) to determine the extent of leukocyte infiltration. Dotted lines demarcate the aortic walls (AW). L—lumen. Scale bar: 20 µm. Right: Box‐whisker plot of mean fluorescence intensity (MFI) of CD45 signals as a readout for the amount of infiltrated CD45+ leukocytes within the aortic walls. Box‐whisker plots indicate median (middle line), 25th, 75th percentile (box), and the lowest/highest data points (whiskers). n = 3 mice per group, with 2 tissue sections analyzed per animal. **P < 0.01 (Mann–Whitney).

Figure EV3

Human plasma and mouse serum studies

Chemokine expressions in the BOECs treated with autologous plasma. Fold changes in gene expressions of CXCL10, CXCL11, CXCL12, CCL20, and CX3CL1 in NAFLD (n = 10) BOECs were normalized to healthy (n = 12) BOECs. Results are indicated by mean ± SD. *P < 0.05; **P < 0.01; ****P < 0.0001; ns—not significant (Mann–Whitney test).

Plasma concentrations of CXCL10 and CXCL12 in healthy and NAFLD subjects. Box‐whisker plots indicate median (middle line), 25th, 75th percentile (box), and the lowest/highest data points (whiskers). ns—not significant (Mann–Whitney test).

Serum CXCL12 concentrations of wild‐type C57BL/6J mice placed on LIDPAD or control diet for 12 and 16 weeks (n = 5 mice), determined by ELISA. Data presented as box‐whisker plots of median (middle line), 25th, 75th percentile (box), and the lowest/highest data points (whiskers). ns—not significant (two‐way ANOVA).

Previous studies have also shown that serum CXCL10 and CXCL12 levels were statistically higher by the severity of fibrosis, steatosis, and inflammation (Chalin et al, 2018). Our ELISA quantification of these chemokines in patient plasma samples, although insignificant, showed a consistent trend with higher CXCL10 and CXCL12 levels in NAFLD plasma than in healthy controls (Fig EV3B). We then tested whether NAFLD endothelial cells could be a contributor to soluble chemokines. Our ELISA data confirmed an intensified secretion of CXCL12 protein in NAFLD BOECs‐conditioned culture media, after correction against nonconditioned media (Fig 3D).

We decided to bring forward CXCL12 as the key candidate to validate in two murine NAFLD models. We first examined liver sections harvested from humanized mice that were reconstituted with human immune components (Her et al, 2020). The high‐fat and high‐calorie (HFHC) diet induce the development of liver steatosis, accompanied by inflammation and fibrosis, which are characteristic features of human advanced NAFLD. Humanized mice on a 20‐week HFHC diet demonstrated higher levels of CXCL12 associated with the liver vasculatures as compared to mice on chow diet, indicated by co‐localization of CXCL12 with CDH5 (VE‐cadherin)‐expressing endothelial cells in the livers (Fig 3E).

In addition, an improved diet‐inducible liver disease mouse model was established by the use of a modified Liver Disease Progression Aggravation Diet (LIDPAD; Low et al manuscript in submission elsewhere). By weeks 8–12, all mice under LIDPAD developed widespread steatosis, intralobular inflammation with ballooning, and infiltration of lymphocytes and macrophages. This model was able to recapitulate the entire spectrum of NAFLD progression (more details in Materials and Methods). As systemic vasculopathies may be associated with NAFLD, we analyzed aortic tissue sections harvested from mice after 12 weeks of LIDPAD or control diets. Histological staining showed that CXCL12 was ubiquitously expressed throughout the aortic intima and media (representative images, Fig 3F). We identified CDH5‐positive endothelial cells that showed an intercellular junctional staining pattern (CDH5 in red, Fig 3F). Quantitatively, CXCL12 expressions within regions of interest strictly containing the CDH5‐positive endothelial lining were found to be significantly higher in LIDPAD aortic sections compared with the chow diet group. Furthermore, we profiled the systemic levels of CXCL12 in serum samples of mice that had been fed on LIDPAD or chow diets for 12 and 16 weeks. However, we did not resolve significant differences of CXCL12 in the mouse sera between LIDPAD and control groups (Fig EV3C). This could be due to soluble CXCL12 in circulation potentially contributed by several sources of cells, including the liver sinusoidal endothelial cells, cholangiocytes, and adipocytes (Bigorgne et al, 2008; Mendt & Cardier, 2012; Kim et al, 2014). We were motivated to find out whether the intensified endothelial‐associated CXCL12 in NAFLD mice could exert a chemotactic effect locally to recruit immune cells. In our previous study, T lymphocyte infiltration into the liver was apparent in mice on HFHC diet (Her et al, 2020). Likewise, we observed increased retention of CD45+ leukocytes in aortic walls of mice fed with NAFLD‐inducing LIDPAD diet as compared to those on chow diet (Fig 3G). Nonalcoholic fatty liver disease endothelial chemokines could play a role in regulating the infiltration of immune cells into vessel walls.

NAFLD endothelial cells may preferentially interact with effector lymphocytes

Interplay of endothelial cells and leukocytes implicates vascular inflammation and atherogenesis (Libby, 2021). We hypothesized that analysis of chemokine receptor expressions in the patient immune populations would allow us to better predict immune interactors of vascular endothelia in NAFLD. Using single‐cell RNA‐sequencing, we performed immunoprofiling on PBMC samples from 2 NAFLD patients (1 male, 1 female) and 4 healthy controls (2 males, 2 females) (Appendix Table S4). In the UMAP visualization, we identified 14 distinct cell type clusters by DICE annotation (Schmiedel et al, 2018; upper panel, Fig 4A). Based on chemokine ligand‐receptor mapping (Hughes & Nibbs, 2018), we analyzed expressions of the counter receptors (lower panel, Fig 4A) to those chemokine candidates enriched in NAFLD endothelial cells. Most of the chemokine receptor genes showed differential expressions across various immune populations, while CXCR4 was more ubiquitously expressed. The chemokine receptor‐expressing populations were more consistently found in T lymphocytes, natural killer (NK) cells, monocytes, and B cells (Fig 4B). Even though there were increased expressions of candidate chemokines in NAFLD endothelial cells, the expressions of their counter receptors, particularly CXCR4 and CX3CR1, appeared lower in NAFLD PBMCs than in healthy controls (Fig 4A). This might be due to the diminished proportions of NK cells and monocytes in NAFLD (Fig 4C), which were among the main chemokine receptor‐expressing populations. Interestingly, the B cell subset in healthy controls seemed to be expressing more CXCR4 compared to that in NAFLD (Fig 4A). B cell‐specific CXCR4 has been shown to be atheroprotective through its regulation of IgM levels (Doring et al, 2020). Whether lower CXCR4 has a casual role in NAFLD‐related vasculopathy requires further interrogation as CXCR4 can be both protective and pathogenic depending on cell‐specific functions and different stages of atherosclerosis development (Murad et al, 2021). To validate CXCR4 and CX3CR1 protein expression in larger sample size, we conducted magnetic‐activated immune cell sorting to derive CD8+ T cells and CD8− cell fractions (largely include monocytes, CD4+ T cells, NK cells, and B cells) from NAFLD (n = 14) and healthy (n = 6) PBMCs, followed by flow cytometry characterization (Fig 4D). About 20–50% of each immune subset expressed CXCR4, but there was no significant difference between NAFLD and healthy groups. Both CD8+ T cells and CD8− lymphocytes subsets had relatively higher levels of CX3CR1+ cells (30–80%), while CX3CR1 expression was lower in monocytes (10–25%). Overall, there were no significant differences in chemokine receptor expressions in the immune cell subsets between NAFLD and healthy groups.

Figure 4

Immunoprofiling of NAFLD patients and healthy subjects by single‐cell transcriptomics

Upper panel: Single‐cell UMAP of sequenced PBMCs from 2 NAFLD patients and 4 healthy subjects; n = 10,786 cells. We identified 15 distinct cell types by DICE annotation. Lower panel: Expression patterns of CXCR3, CXCR4, ACKR3, CX3CR1, and CCR6 across different cell types in NAFLD and healthy PBMC samples.

Chemokine ligand‐receptor mapping to identify potential immune cell interactors based on NAFLD endothelial chemokine signatures.

Graphic shows the percentage composition of identified cell types in NAFLD and healthy PBMCs single‐cell analysis.

Upper panel: Representative flow cytometry dot plots, characterizing the isolated CD8+ T cells and CD8− immune cell fraction (mainly CD4+ T cells, NK cells, and monocytes) from PBMCs of healthy and NAFLD subjects. Lower panel: Flow cytometry quantification of % of positively expressing CXCR4+ and CX3CR1+ cells within each immune subset. Sample sizes are n = 14 NAFLD; n = 6 healthy. Box‐whisker plots indicate median (middle line), 25th, 75th percentile (box), and the lowest/highest data points (whiskers). ns—not significant (t‐test).

Heatmap of gene expressions of endothelial antigens relevant to interaction with T lymphocytes in NAFLD (n = 3) and healthy (n = 3) BOECs, with 3 biological replicates per donor.

Upper panel: Representative flow cytometry histograms of HLA‐A/B/C, HLA‐DM, and IFN‐γR1 expressions in NAFLD (red) and healthy (blue) BOECs. Gray, isotype control. Lower panel: Box‐whisker plots reflect mean fluorescence intensity (a.u—arbituary unit), normalized to isotype control, as a readout for epitope density. Box‐whisker plots indicate median (middle line), 25th, 75th percentile (box), and the lowest/highest data points (whiskers). ns—not significant (t‐test). n = 4 donors; *P < 0.05; ns—not significant (Mann–Whitney test for HLA‐A/B/C and IFN‐γR1; t‐test for HLA‐DM).

We observed an overarching T cell intensification in NAFLD (Fig 4C). To decipher whether NAFLD endothelial cells could potentially interact with T lymphocytes, we retrieved relevant gene expressions from the BOEC transcriptomic dataset, particularly the human leukocyte antigen (HLA) family members. HLA class I molecules mainly support interaction with CD8+ T cells, whereby HLA class II molecules facilitate binding with CD4+ T subsets (Laidlaw et al, 2016). We found significantly higher expressions of several HLAs in NAFLD BOECs than in healthy BOECs, involving class I molecules HLA‐A, HLA‐E, and HLA‐H and class II molecules HLA‐DMA and HLA‐DMB (Fig 4E). NAFLD BOECs also had significantly higher expressions of interferon‐γ receptors IFNGR1 and IFNGR2 (Fig 4E), potentially making them susceptible to interferon‐γ that can modify the activated phenotype of endothelial cells to favor T helper 1 cell inflammatory reactions (Pober & Sessa, 2007). We then performed flow cytometry to evaluate these markers at the protein levels. Almost all NAFLD and healthy BOECs expressed HLA class I antigens (HLA‐A/B/C), but there were minimal expressions for HLA‐DM class II antigens (Fig 4F, upper panel). We did not notice major differences in IFN‐γR1 expressions between NAFLD and healthy BOECs. Interestingly, epitope density of HLA class I antigens was significantly higher in NAFLD BOECs than in healthy controls (Fig 4F, lower panel: mean fluorescence intensity). The data thus far suggest that endothelial cells in NAFLD may be providing a chemokine milieu and partially enhanced antigen presentation, which could potentially promote their interaction with immune cells.

The CXCL12‐CXCR4 axis is implicated in NAFLD endothelial‐immune crosstalk

Phenotypically, to study the interaction of NAFLD endothelial cells with immune mediators, we performed coculture assays of endothelial‐immune cell adhesion, chemotaxis, and transendothelial migration. We first employed two types of human leukocyte cell lines—T lymphocytic Jurkat cells and monocytic THP‐1 cells for adhesion assay. Jurkat cells demonstrated a greater extent of adhesion onto both NAFLD and healthy BOEC monolayers than THP‐1 cells (Fig 5A). Significantly more Jurkat cells also adhered to NAFLD BOECs than to healthy BOECs. Since NAFLD BOECs were found to be activated with chemokine hallmarks (Fig 3), we further confirmed that NAFLD BOECs had significantly greater chemotactic ability than healthy BOECs in recruiting Jurkat cells but not THP‐1 cells (Fig 5B). These data seemed to indicate preferential recruitment of T lymphocytes to NAFLD BOECs, compared with monocytes.

Figure 5

Characterizations of endothelial‐immune interactions

Left: Workflow schematic of adhesion assay. Middle: Representative images of labeled THP‐1 and Jurkat cells attached to BOEC monolayers (Scale bars, 100 µm). Right: Box‐whisker plots of the numbers of adhered leukocytes were determined fluorometrically. Sample sizes are n = 9 NAFLD and n = 7 healthy, with 3 biological replicates per donor. Box‐whisker plots indicate median (middle line), 25th, 75th percentile (box), and the lowest/highest data points (whiskers). *P < 0.05; ****P < 0.0001; ns—not significant (two‐way ANOVA).

Left: Schematic of chemotaxis assay. Right: Box‐whisker plots of the numbers of recruited THP‐1 and Jurkat cells by the BOECs. Box‐whisker plots indicate median (middle line), 25th, 75th percentile (box), and the lowest/highest data points (whiskers). ns—not significant (t‐test). Sample size is n = 5 donors; *P < 0.05; ns—not significant (t‐test).

CD8+ T cells and CD8− cell fractions were isolated from NAFLD and healthy subjects. Box‐whisker plots show the number of recruited immune cells to their respective BOECs matched by health condition in coculture chemotaxis assays. Box‐whisker plots indicate median (middle line), 25th, 75th percentile (box), and the lowest/highest data points (whiskers). ns—not significant (t‐test). Sample sizes are n = 11 NAFLD and n = 6 healthy, with 3–4 technical replicates per donor. *P < 0.05; ns—not significant (Mann–Whitney).

Box‐whisker plots show percentages of inhibition on immune cells chemotaxis by treatment with AMD3100 (CXCR4 small molecule inhibitor). Box‐whisker plots indicate median (middle line), 25th, 75th percentile (box), and the lowest/highest data points (whiskers). ns—not significant (t‐test). Sample sizes are n = 6 NAFLD and n = 4 healthy, with 3–4 technical replicates per donor. ****P < 0.0001; ns—not significant (two‐way ANOVA).

Box‐whisker plots show percentages of inhibition on immune cells chemotaxis by treatment with anti‐CXCL12 IgG (CXCL12‐neutralizing antibody). Box‐whisker plots indicate median (middle line), 25th, 75th percentile (box), and the lowest/ highest data points (whiskers). ns—not significant (t‐test). Sample sizes are n = 5 NAFLD and n = 3 healthy, with 3–4 technical replicates per donor. ****P < 0.0001; ns—not significant (two‐way ANOVA).

Left: Schematic of transendothelial migration assay with patient‐derived BOECs and immune cells. Right: Transendothelial electric resistance (TEER) measurements of BOEC barrier tightness in the presence and absence of PBMCs and/or AMD3100. Box‐whisker plots indicate median (middle line), 25th, 75th percentile (box), and the lowest/highest data points (whiskers). ns—not significant (t‐test). Sample size is n = 5 patient samples, with 3 technical TEER readings recorded per sample. ****P < 0.0001; ns—not significant (one‐way ANOVA).

Box‐whisker plots of the numbers of circulating endothelial cells per million PBMCs in healthy (n = 6) and NAFLD subjects (n = 11). Box‐whisker plots indicate median (middle line), 25th, 75th percentile (box), and the lowest/highest data points (whiskers). ns—not significant (t‐test). *P < 0.05 (Mann–Whitney test).

To better mimic the pathophysiological context, we switched to using PBMCs obtained from NAFLD patients and healthy subjects for further co‐culture studies with BOECs. We conducted magnetic‐activated immune cell sorting to derive CD8+ T cells and CD8− cell fraction (largely including CD4+ T cells, NK cells, monocytes) from NAFLD and healthy PBMCs. Functionally, we assessed the chemotactic migration of NAFLD and healthy CD8+/− immune cell fractions toward BOEC monolayers matched by their health conditions. Compared with healthy BOECs, NAFLD BOECs were able to recruit significantly more CD8+ T cells, while no apparent difference was observed for the chemotaxis of CD8− cell fraction (Fig 5C). This could be in consensus with NAFLD BOECs possessing higher epitope density of HLA class I antigens (Fig 4F) that mediate interaction with CD8+ T cells.

Since CXCL12 was consistently enhanced in NAFLD patient endothelial cells, and in vascular endothelia of in vivo disease models (Fig 3), we postulated the involvement of the CXCL12‐CXCR4 axis in the immune cell infiltration of blood vessels in NAFLD. To impair the CXCL12‐CXCR4 axis in endothelial‐immune chemotaxis, we first tested AMD3100, a small molecule inhibitor of CXCR4, as it has been validated and approved for human use (Hatse et al, 2002; Fricker et al, 2006). The specificity of AMD3100 is attributed to its well‐defined high‐affinity molecular interactions between the cyclam moieties of drug molecules and aspartate residues of the CXCR4 receptor (Gerlach et al, 2001). Here, we found that chemotactic migration of CD8− cell fraction was attenuated more significantly than CD8+ T cells in the presence of AMD3100 in NAFLD coculture assay but not in the healthy group (Fig 5D). This was despite no major difference in CXCR4 expressions on the various immune cell subsets of NAFLD and healthy subjects, which we described in Fig 4D. Besides chemokine receptor expression, it was still unknown to us whether the chemokine binding affinity of the CXCR4 receptor could be different in various immune subsets or altered under disease context as reported in other studies (Bigorgne et al, 2008; Boujedidi et al, 2015). Next, we tested a CXCL12‐neutralizing antibody to counteract the effect of enhanced chemokine secretion from NAFLD BOECs. Indeed, the use of anti‐CXCL12 IgG led to a significantly greater inhibition of NAFLD CD8+ T cell chemotaxis than in the healthy group (Fig 5E). Anti‐CXCL12 IgG was similarly effective at moderating the chemotaxis of CD8− cell fraction in both NAFLD and healthy groups. From Fig 5D and E, we were mindful about a small extent of increased migration of immune cells (shown as negative data points of % inhibition) in the presence of AMD3100 or anti‐CXCL12 IgG. It has been described that inhibition of the CXCL12‐CXCR4 axis may cause a shift in cytokine profiles of T cells (Chen et al, 2019; Zeng et al, 2019). We postulate that any potential change of cytokine profiles in T cells because of AMD3100 or anti‐CXCL12 IgG treatment may in turn affect our BOEC activation and dynamics of immune cell chemotaxis. While the CXCL12‐CXCR4 axis may be implicated in NAFLD endothelial‐immune crosstalk, further studies in preclinical models are necessary to establish the effectiveness of blocking either the chemokine or its chemokine receptor.

To understand the vascular consequences of endothelial‐immune crosstalk, we characterized NAFLD BOEC barrier integrity after transendothelial migration assay with NAFLD PBMCs. We measured transendothelial electrical resistance (TEER) of BOEC monolayers, where lower TEER readings would indicate greater endothelial barrier permeability. At baseline (Fig 5F, Untreated group), we did not observe changes in TEER upon transmigration of NAFLD PBMCs. However, perturbations of the endothelial barrier became apparent when NAFLD BOECs were pretreated with LPS and FFA, as seen in the significant reduction of TEER readings after interaction with NAFLD PBMCs (Fig 5F, LPS + FFA group). Both LPS and FFA are factors known to be implicated in NAFLD, potentially rendering BOECs more sensitive to immune‐inflicted injury. Inhibition of the CXCL12‐CXCR4 axis by AMD3100 averted the impact of PBMC transmigration and restored TEER readings (Fig 5F, LPS + FFA group). Taking these findings to the clinical setting, we analyzed circulating endothelial cells (CECs) that were cells shed from damaged blood vessels into the bloodstream, hence indicative of the level of vascular injury (Hebbel, 2017). NAFLD patients had more than 3 folds of CECs compared with healthy subjects (Fig 5G), validating that endothelial pathology was indeed aggravated under NAFLD condition.

Discussion

With the prevalence of NAFLD, there is an increasing recognition of cardiovascular complications being significant causes of mortality in these patients. In this study, we have advanced mechanistic understanding of endothelial pathophysiology in NAFLD, involving an interplay with immune cells (Fig 6). Our key findings are: 1. NAFLD endothelial cells were activated and enriched for chemokine hallmarks and probably some human leukocyte antigens; 2. the CXCL12‐CXCR4 axis may be preferentially implicated in the interaction of NAFLD endothelial cells with immune cells; and 3. NAFLD patients had a more pronounced endothelial injury, which may in part be due to heightened interactions with effector immune cells.

Figure 6

Mechanisms of vascular injury due to NAFLD endothelial chemokine milieu and preferential interaction with effector immune cells via the CXCL12‐CXCR4 axis.

Several studies have reported elevated chemokines in the plasma/sera of NAFLD patients, and the role of these chemokines in regulating liver inflammation (Marra & Tacke, 2014; Pan et al, 2020). Potential sources of cells capable of contributing to soluble CXCL12 in circulation include liver sinusoidal endothelial cells, cholangiocytes, and adipocytes (Bigorgne et al, 2008; Mendt & Cardier, 2012; Kim et al, 2014). Here, the transcriptomic analysis demonstrated enriched chemokine signatures in NAFLD patient endothelial cells, suggesting that activated endothelial cells could be a potential source of elevated chemokines in NAFLD patients. This was supported by increased in vivo expressions of CXCL12 associated with arterial endothelial cells and liver vasculatures in our mouse NAFLD models. Although CXCL12 has previously been studied for its beneficial impact on atherosclerotic plaque stability and neurovascular repair after stroke, these mainly relate to cell‐specific atheroprotective effects and CXCL12’s role in the recruitment of bone marrow‐derived cells to sites of injury (Wang et al, 2012; Zernecke & Weber, 2014). Instead, plasma CXCL12 levels predict long‐term outcomes for stroke and adverse cardiovascular events (Schutt et al, 2012; Ghasemzadeh et al, 2015). Mechanistically, CXCL12 is known to disrupt cholesterol efflux from macrophages by activating GSK‐3β/β‐catenin/TCF21 signaling pathway, thus exacerbating atherosclerosis (Gao et al, 2019). Notably, endothelial‐derived CXCL12 is a key driver of atherosclerosis and one of the contributors to plasma CXCL12 levels (Döring et al, 2019). Atherosclerosis Apoe −/− mice deficient of CXCL12 from arterial endothelial cells, but not from smooth muscle cells, bone marrow–derived hematopoietic, and nonhematopoietic resident cells, led to a marked reduction of lesion size (Döring et al, 2019). Given that chemokine secretions by endothelial cells could influence the recruitment of distinct immune populations (Øynebråten et al, 2004), we further postulated that NAFLD endothelial cell‐derived chemokines could act in a paracrine manner to recruit immune cells, possibly leading to immune‐mediated endothelial damage.

While the immune landscape of NAFLD is not yet a well‐studied area, there are indications of immune cell population changes and their contribution to NAFLD progression. Alterations of cytotoxic CD8+ T cells are associated with steatohepatitis phenotypes in patients, and CD4+ T helper cells (Th1 and Th17) are positively correlated with the severity of NAFLD (Haas et al, 2019). In our previous work using humanized mice reconstituted with human immune components, we found that CD4+ T cells promoted liver fibrosis (Her et al, 2020). As obesity is often associated with NAFLD, obesity‐induced adipose tissue inflammation has been shown to involve CXCL12 secreted from adipose tissue, which recruits macrophages and induces insulin resistance in mice fed with a high‐fat diet (Kim et al, 2014). Here, our single‐cell analysis pointed to an overarching T cell intensification in NAFLD patients. Correspondingly, our coculture assays using patient‐derived endothelial and immune cells revealed preferential chemotaxis of the CD8+ T cell subset, and increased endothelial permeability after immune cell transmigration. Endothelial cells are known to be one of the primary targets of immunologic attack. Cytotoxic CD8+ T cells have been shown to promote endothelial cell apoptosis in the vulnerable atherosclerotic plaques of Apoe −/− mice (Kyaw et al, 2013). In patients with acute coronary syndrome, cytotoxic effector molecules from CD8+ and CD4+ T cells cause damage to endothelial cells and hence plaque destabilization (Nakajima et al, 2002; Leistner et al, 2020). Further studies will be needed to elucidate how T cell subsets regulate inflammation in the NAFLD setting.

We believe that activated endothelial cells in NAFLD may be providing a chemokine milieu that preferentially interplays with effector immune cells through chemokine ligand‐receptor interactions, potentially leading to inflammatory pathology and cytotoxicity‐induced vascular injury. The CXCL12‐CXCR4 axis has been shown to mediate intrahepatic infiltration of immune cells in mouse diet‐induced models of obesity and NASH, where it has been postulated that change in membrane cholesterol may alter CXCR4 receptor conformation, hence increasing its chemokine binding affinity (Bigorgne et al, 2008; Boujedidi et al, 2015). These studies report no change in either liver expression of CXCL12 or intrahepatic T cell expression of CXCR4. Also, it has been shown in their studies that liver CXCL12 originates mainly from bile duct epithelial sources and not from endothelial cells. Here, our work demonstrated an endothelial source of CXCL12 in patient cells, and in the aortae and liver vasculatures of NAFLD mouse models. Since we did not notice significant differences in CXCR4 expressions in NAFLD and healthy immune cells, the upregulation of endothelial CXCL12 in NAFLD might be primarily responsible for immune cell recruitment to vasculatures. We recognize the limitation that we had not investigated the chemokine binding affinity of CXCR4 receptors on our patient immune cells.

Current therapeutic exploration on the intervention of chemokines and their receptors in NAFLD are mainly targeting hepatocellular carcinoma (Chen et al, 2018). The recent outcomes of the Canakinumab Anti‐inflammatory Thrombosis Outcome Study (CANTOS) prove that immunomodulatory treatment can be used to treat atherosclerosis by reducing the inflammatory burden. However, such treatment that uses monoclonal antibody targeting interleukin‐1β is impeded by its high price, failure to improve overall mortality, and an increased risk of infection (Ridker et al, 2017; Sehested et al, 2019). Furthermore, the limited success of the Cardiovascular Inflammation Reduction Trial (CIRT) that tested low‐dose methotrexate (Ridker et al, 2019), has shown that broad anti‐inflammatory treatments are ineffective in reducing cardiovascular events. Therefore, new affordable therapeutics tailored to specific autoimmune mechanisms and with limited immunosuppressive effects are still required. Currently, the CXCR4 inhibitor, AMD3100, is used clinically for the treatment of HIV infection and as a hematopoietic stem cell (HSC) mobilizer in non‐Hodgkins lymphoma (NHL) patients receiving radiation therapy (De Clercq, 2005; Keating, 2011). More recently, experimental therapeutics using AMD3100 have been explored in Alzheimer’s disease (AD; Gavriel et al, 2020), amyotrophic lateral syndrome (ALS; Rabinovich‐Nikitin et al, 2016), and hepatopulmonary syndrome (HPS; Shen et al, 2018). In these conditions, AMD3100’s interference with the CXCR4‐CXCL12 signaling axis proves useful in reducing inflammatory burden, resulting in attenuated disease progression in vivo. Nonetheless, CXCR4 can be both protective and pathogenic depending on cell‐specific functions and different stages of atherosclerosis development (Doring et al, 2020; Murad et al, 2021). To counteract the chemotactic effect of NAFLD endothelial‐derived chemokine, we found that treatment with a CXCL12‐neutralizing antibody significantly reduced NAFLD CD8+ T cell chemotaxis compared to that of the healthy group in our endothelial‐immune coculture assays. It has been found that CXCL12 inhibition by chalcone 4, a small neutralizing compound binding to CXCL12, confers greater efficacy than CXCR4 blockade with AMD3100 in decreasing macrophage infiltration in the lungs of rat models of pulmonary hypertension (Bordenave et al, 2020). Another RNA oligonucleotide that neutralizes CXCL12, known as NOX‐A12, is currently under development in multiple oncology indications. NOX‐A12 has been found to interfere with chronic lymphocytic leukemia cell migration and drug resistance (Hoellenriegel et al, 2014). While drug repurposing to impair the CXCL12‐CXCR4 axis may be a viable avenue to mitigate NAFLD‐associated vasculopathies, it is important for further preclinical testing to establish the effectiveness of blocking chemokine and/or its receptor.

Importantly, we found that inhibition of the CXCL12‐CXCR4 axis could reduce endothelial barrier impairment posed by FFAs/LPS exposure and transendothelial migration of immune cells. In evaluating vascular injury in NAFLD, this is the first report of a significantly higher level of circulating endothelial cells (CECs) in NAFLD patients than in healthy controls. CECs are not to be confused with endothelial progenitor cells (EPCs), which originate from bone marrow and are recruited to sites of vascular injury to undergo repair (Hebbel, 2017). An altered status of circulating EPCs hence reflects one’s endothelial repair capacity. NAFLD patients were found to have reduced levels of circulating EPCs and impaired adhesive and migratory functions than their healthy counterparts. This suggests that the attenuation of EPC‐mediated endothelial repair may contribute to atherosclerotic disease progression and enhanced cardiovascular risk and events in NAFLD patients (Chiang et al, 2012; Francque et al, 2016). Circulating endothelial cells, on the other hand, are directly indicative of in situ vascular injuries, as they are cells dislodged from damaged endothelial into the bloodstream (Hebbel, 2017). A limitation of this study is that we have not distinguished NAFLD patients with and without confounding cardiovascular risks such as hypertension and diabetes, which may also contribute to more CECs (Chioh et al, 2021). CECs are a representative marker of endothelial dysfunction and are useful to serve as a form of vascular health surveillance in conditions that predispose patients to vascular complications.

We present that vasculopathy in NAFLD may involve endothelial chemokine activation and vascular‐immune interactions. Vascular inflammation may be further exacerbated by shared cardiometabolic risk profiles common to most NAFLD patients. While there is much complexity in the pathogenesis of NAFLD, impairing specific pathways of endothelial‐immune crosstalk may be effective at reducing the burdens of vasculopathy in NAFLD patients. Finally, it is imperative to monitor the vascular health risk of thromboembolism as part of the management of NAFLD patients so as to ensure early preventive therapies.

Materials and Methods

Study approvals, patient enrollment, and sample collection

This study was approved by the Local Ethics Committee of National Healthcare Group Domain Specific Review Board (DSRB Ref: 2016/00580) and Nanyang Technological University Singapore Institutional Review Board (IRB‐2020‐05‐012 and IRB‐2020‐09‐011). Written informed consent was obtained from each participant after the nature, and possible consequences of the studies have been explained. The study protocol complies with the Helsinki Declaration.

Patients with biopsy‐proven NAFLD were taken from the fatty liver clinic or from the bariatric surgery clinic. Patients with NAFLD had at least 5% steatosis characterized by liver biopsies. NAFLD patients were selected across the spectrum of disease. Controls were selected from healthy volunteers having a Controlled Attenuation Parameter score of less than 248 on vibration controlled transient elastography, indicative of less than 5% steatosis. In addition, control patients did not have raised liver enzymes or evidence of fat in the liver on any other imaging modality.

For sample collection, 10 ml of fresh blood was collected from each participant and processed in the laboratory within 6 h. Upon Ficoll centrifugation of the blood specimen, a buffy coat layer containing peripheral blood mononuclear cells (PBMCs) was isolated. PBMC fractions were used for three purposes—(i) cultivated in cell culture to derive blood outgrowth endothelial cells (BOECs); (ii) immunophenotyping and isolation of immune subpopulations; and (iii) analysis of circulating endothelial cells.

Derivation and maintenance of blood outgrowth endothelial cells

We employed established protocols to derive our BOECs (Martin‐Ramirez et al, 2012a; Ormiston et al, 2015). Peripheral blood samples (5–9 ml) were diluted 1:1 with phosphate‐buffered saline and separated to obtain a buffy coat by density gradient centrifugation over Ficoll® Paque (GE Healthcare). The buffy coat, which was enriched with PBMCs, was carefully collected, washed with PBS, resuspended in endothelial growth media 2 (EGM‐2) medium (Lonza) supplemented with 16% defined fetal bovine serum (FBS; Hyclone). We typically seeded 20–22 million PBMCs from each donor onto collagen I‐coated wells at a cell density of 2 × 106 cells/cm2. The culture medium was changed every two to three days. Outgrowth colonies should then appear between 7‐ to 14‐day postseeding. We obtained an average of 1–2 blood outgrowth endothelial colonies per 10 million cultivated PBMCs. Quality controls on cell purity, endothelial marker expressions, and functional characterization were performed on BOECs from passages 2 and 3. Only BOEC lines that passed quality control checks would be utilized for downstream experimentations.

Endothelial marker characterization by flow cytometry

Cell surface markers were quantified using flow cytometry to phenotypically confirm the endothelial identity of the derived BOECs. Briefly, BOECs were washed and stained with a cocktail of surface antibodies in a staining buffer containing PBS with 2% FBS, for 15 min in the dark. CD31 (platelet endothelial cell adhesion molecule 1; clone wm59, BioLegend) and CD144 (vascular endothelial cadherin; clone 55‐7H1, BD Pharmingen) were selected as endothelial markers CD45 (clone 2D1, Invitrogen) as leukocyte exclusion marker and CD133 (clone 7, BioLegend) as a progenitor marker. Fluorescence data were collected on a BD LSR Fortessa X‐20 cell analyzer (Becton Dickinson) and analyzed using FlowJo v10.7.1 software (Becton Dickinson). Details on antibodies are included in Appendix Table S5.

Endothelial marker characterization by immunocytochemistry

Confluent monolayers of BOECs were fixed and stained with anti‐CD31 (clone 89C2, 1:200, Cell Signaling) and anti‐CD144 (clone C‐19, 1:200, Santa Cruz) for identity verification. Staining was detected with secondary antibodies Alexa‐Fluor 488 and 568 and visualized using ZEISS Celldiscoverer7 microscope system. All experiments were accompanied by a negative control—fibroblasts. Details on antibodies are included in Appendix Table S5.

Endothelial functional characterization by angiogenesis assays

Fibrin gel bead‐based sprouting angiogenesis assay was performed as per described (Nakatsu et al, 2007). Briefly, BOECs were coated onto collagen‐coated Cytodex 3 microcarrier beads (Sigma‐Aldrich) at 150 cells/bead for 4 h and allowed to adhere overnight. The coated beads were then suspended in a 2.0 mg/ml of fibrinogen solution (Sigma‐Aldrich) at a concentration of 500 beads/ml. Fibrin gels were formed by adding 500 µl of the fibrinogen/bead suspension to each well of a 24‐well plate containing 0.625 U/ml of thrombin (Sigma‐Aldrich). Once clotted, the gels were topped with 0.5 ml of EMG‐2 supplemented with 10% heat‐inactivated FBS (Gibco). Cells were incubated at 37°C/5% CO2, and sprouting was observed after 24 h using bright‐field microscopy (Lumenera).

BOEC RNA‐sequencing and analysis

Total RNA was isolated from the BOECs as per the manufacturer’s instructions (RNeasy Micro Kit, Qiagen). PolyA library preparation and 150 base‐pair paired‐end sequencing was performed by the Novogene sequencing facility (Singapore) using an Illumina HiSeq sequencer. The average sequencing depth was 50 million reads. All subsequent analyses were performed on Partek Flow software, version 7.0 (St. Louis, MO, USA). Reads were aligned to the human reference genome GRCh37 using STAR version 2.6.1d. Transcript abundances were quantified with the Partek E/M algorithm based on Ensembl Transcripts Release 83, subsequently normalized by CPM (counts per million) method. Differentially expressed genes in the transcriptome data were identified using the gene‐specific analysis (GSA) method for P‐values and the Benjamini–Hochberg method for FDR. Genes were called differentially expressed whether P‐value < 0.05, FDR < 0.05, and fold change ≥ 2 or ≤ −2. Enrichment analyses for differentially expressed genes were performed using Metascape (Zhou et al, 2019). Protein‐protein interactions were analyzed using the STRING tool, and the yielded networks were further constructed and analyzed in Cytoscape (version 3.8.0) with the network topography based on the confidence score calculated by STRING. Cytoscape plugin cytoHubba was used to assess the centrality (strongly interconnected subnetworks) of differentially expressed genes in the protein‐protein interaction networks with 12 algorithms: Maximum Neighborhood Component (MCC), Density of Maximum Neighborhood Component (DMNC), Maximal Clique Centrality (MNC), Degree (Deg), Edge Percolated Component (EPC), BottleNeck (BN), EcCentricity (EC), Closeness (CLO), Radiality (Rad), Betweenness (BC), Stress (STR), and ClusteringCoefficient (CC). The selected hub genes were verified by the MCODE method with default algorithms (degree cut‐off of 2, node score cut‐off of 0.2, K‐Core of 2, and maximum depth of 100).

Quantitative PCR

Total RNA was isolated using RNeasy Plus Mini kit (Qiagen) as per the manufacturer’s protocol, subsequently used to generate cDNA with LunaScriptTM RT SuperMix Kit (New England Biolabs). Real‐time PCR was performed using SYBR green gene expression assays (New England Biolabs) on a QuantStudio 6 instrument (Applied Biosystems). Gene expressions were normalized to the endogenous GAPDH housekeeping gene.

Flow cytometry

Flow Cytometry was performed using BD LSRFortessa X‐20 (BD Biosciences), and data were analyzed using FlowJo v10.7.1 software (Becton Dickinson). Cell counts and Mean Fluorescence Intensity (MFI), indicative of target protein‐expressing cells and ligand abundance, respectively, were computed and normalized to unstained controls prior to statistical comparisons between both study groups. (i) Characterization of chemokine receptors on immune cells. Thawed PBMCs were enriched for CD8 using magnetic‐associated cell sorting (MACS) and subjected to flow cytometry to determine the expression of immune cell markers CD4, CD8, and chemokine receptors CXCR4, CX3CR1. Sorted immune cells were fixed before staining fluorophore‐conjugated antibodies (Appendix Table S5) specific to the abovementioned immune cell markers. Flow cytometry was performed, and data were analyzed using FlowJo as per described above. (ii) Determining surface antigens on BOECs. Confluent monolayers of BOECs collected after 48 h were subjected to flow cytometry to determine the expression of target proteins HLA‐A/B/C, HLA‐DM, and IFNγR. BOECs were fixed and permeabilized before staining with fluorophore‐conjugated antibodies (Appendix Table S5) specific to the abovementioned targets. Flow cytometry was performed, and data were analyzed using FlowJo as per described above.

Treatment of BOECs

Blood outgrowth endothelial cells were grown to confluence in EGM‐2 media supplemented with 10% of their autologous plasma or NAFLD‐related risk factors (100 ng/ml of lipopolysaccharide and/or free fatty acids including 100 µM of palmitic acid and 100 µM of oleic acid) for 48 h. At the end of treatment, cells were lysed for RNA extraction and qPCR. ELISA of chemokines in plasma and BOEC‐conditioned media are described below.

ELISA of chemokines in BOEC‐conditioned media

EGM‐2 media conditioned by confluent monolayers of BOECs were collected after 48 h. Secretion of IP‐10 (CXCL10, Abcam, ab173194), SDF‐1 (CXCL12, Abcam, ab100637), MIP‐3 alpha (CCL20, Abcam, ab269562), and fractalkine (CX3CL1, Abcam, ab192145) into conditioned EGM‐2 were measured by ELISA as per the manufacturer’s protocols. Reads were corrected for background with nonconditioned EGM‐2.

Animal studies approval

Animals' care was in accordance with institutional guidelines, approved by the local Institutional Animal Care and Use Committee (IACUC #: 181367, A18031, A18033). This research complied with the Guidelines on the Care and Use of Animals for Scientific Purposes of the National Advisory Committee for Laboratory Animal Research of Singapore (NACLAR) and the US National Institute of Health (NIH).

Diet‐induced NAFLD murine models

Wild‐type C57BL/6J male mice (8–9 weeks old) were placed on LIDPAD diet (n = 3) or control diet (purified diet containing 14% protein, 76% carbohydrates, 10% fat) for 12 weeks. The manuscript (Low et al, data not published) relating to this novel LIDPAD diet‐induced NAFLD model is in submission elsewhere. This is a diet‐based model, in which NAFLD phenotype is induced by a Liver Disease Progression Aggravation Diet (LIDPAD; Teklad Custom Diet, TD.140897, Envigo, USA) that contains 43.3 kcal% of fat. Briefly, wild‐type C57BL/6J mice (Invivos, Singapore) aged between 8–9 weeks, were given ad libitum to either LIDPAD or control diet (Teklad Custom Diet, TD. 140806, Envigo, USA; 10.3 kcal% of fat) for 12 weeks. The mice displayed a sustained weight gain in both body and liver as early as one week into the LIDPAD diet. This was also accompanied by impaired glucose metabolism, which developed into glucose intolerance within 8 weeks and insulin resistance within 12 weeks. Histologically, steatosis was evident in the livers from weeks 1–4 of LIDPAD. By weeks 8–12, all mice under LIDPAD developed intralobular inflammation with ballooning and infiltration of lymphocytes and macrophages, in addition to more widespread steatosis. Fibrosis was also observed in 80% of the mice by week 12. Together, this model recapitulates the entire spectrum of NAFLD progression including leukocyte infiltration. We did not use any anesthetic agents during the course of the study or endpoint. Mice were euthanized using CO2. Descending thoracic aorta was then harvested and aortic specimens were fixed overnight in 4% paraformaldehyde (PFA) and embedded in paraffin. Every 1st, 10th, and 20th serial sections (5 μm thickness) were mounted onto each slide for immunostaining experiments.

Humanized NAFLD murine models

Humanized mice containing the human immune system (HIL, age 32 weeks) were obtained from the Institute of Molecular and Cell Biology (IMCB), Agency for Science, Technology and Research, Humanized Mouse Unit. Briefly, human immune reconstitution of NOD‐scid Il2rγ (NSG) mice was performed via intrahepatic injection of human CD34+ fetal liver cells into one to three‐day‐old NSG pups following sub‐lethal irradiation. Murine models of NAFLD were generated from mice with more than 10% human immune reconstitution, determined according to described criteria (Her et al, 2020). We did not use any anesthetic agents during the course of the study or endpoint. HIL mice were euthanized using CO2 after being placed on HFHC diet (Surwit high‐fat diet; 58 kcal% fat that is mainly saturated, with carbohydrate‐enriched drinking water) or chow diet for 20 weeks. After 20 weeks of diet, HFHC‐HIL mice demonstrated key pathologies representative of human NAFL and NASH, along with visibly increased leukocyte infiltration relative to chow diet‐fed mice. Serial liver sections (5 μm thickness) harvested from 4 HFHC‐HIL mice (1 male, 3 female) and 4 chow diet‐fed mice (2 male, 2 female) were used for immunostaining experiments.

Immunohistochemical staining and image analysis

Double immunostaining to investigate (i) CXCL12 (SDF‐1) expressions in the endothelial cells (by CDH5 and SDF‐1) or (ii) leukocyte infiltration (by CD45 and CD3) was performed on the tissue sections from both murine models. Briefly, heat‐induced epitope retrieval (HIER) was performed on deparaffinized sections in the R‐Universal epitope recovery buffer (Aptum Biologics, #AP0530‐500), using the 2100 Antigen Retriever (Aptum Biologics). Sections were blocked with 1X PBS with 1% BSA for 1 h to minimize nonspecific staining before being incubated in a primary antibody mixture of either a goat polyclonal anti‐CDH5 (Santa Cruz Biotechnology, sc‐6458; 4 µg/ml) or a rabbit polyclonal anti‐SDF‐1 (Bioss, bs‐4938R; 5 µg/ml); or a rat monoclonal anti‐CD45 (Tonbo, 70‐0452; 5 µg/ml) or a rabbit polyclonal anti‐CD3 (Invitrogen, PA5‐102404; 10 µg/ml) for 1 h at room temperature. The antibody combinations were then detected by a mixture of conjugated donkey secondary antibodies: an AF568 anti‐goat (Invitrogen, A11057; 4 µg/ml) and an AF488 anti‐rabbit (abcam, ab150073; 10 µg/ml), or an DyLight 650 anti‐rat (Invitrogen, SA5‐10029; 2 µg/ml) and an AF568 anti‐rabbit (Invitrogen, A10042; 4 µg/ml), respectively, for 1 h at room temperature. All sections were counter‐stained with DAPI. Slides were mounted using Fluoromount™ Aqueous Mounting Medium (Sigma‐Aldrich, F‐4680) in preparation for image acquisition. All steps were replicated twice and accompanied by a negative control (no primary antibodies) and a splenic positive control. Details on antibodies are included in Appendix Table S5.

For analysis of endothelial‐associated CXCL12, images were acquired at 40× magnification with a ZEISS Celldiscoverer 7 microscope system. Analysis was carried out with Fiji software by two independent researchers (one blinded), on the delimitated vascular endothelia of murine aortae and liver sections (Schindelin et al, 2012). For images of the murine aorta, vascular endothelia were manually delimitated using region of interest (ROI) drawing tools based on CDH5 signals and internal elastin lamellae. Following background subtraction, CXCL12 fluorescent intensity within the delimitated endothelia was quantified as the sum of pixel intensities in the CXCL12 single‐channel image, normalized by the number of pixels with above‐zero intensities in the DAPI single‐channel image. For images of the murine liver, sinusoidal endothelial regions were delimitated by applying histogram‐based automated thresholding on the CDH5 single‐channel image. Using the EzColocalization plugin on Fiji, Mander’s Coefficient (M2) was obtained as a measure for the degree of CXCL12 overlap with CDH5 within delimitated endothelia (Stauffer et al, 2018). M2 indicates the total number of above‐threshold CXCL12 pixels as a fraction of above‐threshold CDH5 pixels. Threshold values (FT) were set to 0.1 for both CXCL12 and CDH5 channels, meaning that only the 10% of pixels with the highest signals in each channel would be analyzed. For both vascular endothelia of murine aortae and liver sections, there were 10–20 independent regions of interest analyzed per mouse for image quantification.

For aortic leukocyte infiltration, analysis was performed on Fiji. Images of complete aortic ring section (magnification of 5×) were divided into four even areas for quantification of CD45+ cells. Aortic walls (including media and adventitia) were delimitated and the signal for CD45 within the regions of interest were then measured as mean fluorescence intensity (MFI). All readings were subtracted for background MFI from a nonfluorescent area of the same image. There were 8–10 independent regions of interest analyzed per mouse.

Immunoprofiling by single‐cell RNA‐sequencing

Cryopreserved PBMCs were thawed in a 37°C water bath and resuspended in Iscove’s modified Dulbecco’s medium (Gibco) supplemented with 10% heat‐inactivated FBS. PBMCs were then pelleted out of the solution by centrifuging at 300 g for 5 min. Cell number and viability percentage were determined using Trypan blue and an automated cell counter (Countess II, Thermo Fisher Scientific). NAFLD PBMCs were pooled from 2 NAFLD patients while healthy PBMCs were pooled from 4 healthy individuals with each individual contributing about 200,000 cells each. Final PBMCs pools were pelleted down at 300 g for 5 min and resuspended in PBS with 0.04% BSA to give a final concentration of ~1,200 cells/μl.

Cell suspensions were then loaded onto a 10× Genomics Chromium Controller chip by facility personnel at Single‐cell Omics Centre (SCOC), Genome Institute Singapore (GIS). NAFLD and healthy PBMC pools were prepared as separate scRNA‐seq libraries using Chromium Single Cell 3’ v3 Reagent Kit (10× Genomics) by SCOC GIS, and the final ready‐to‐sequence libraries were handed over with quantification and quality assessment reports from Bioanalyzer Agilent 2100 using the High Sensitivity DNA chip (Agilent Genomics). NAFLD and healthy PBMC libraries were pooled equimolarly and sent for sequencing by NovogeneAIT Genomics (Singapore). Raw sequencing data were also processed by NovogeneAIT Genomics (Singapore) using CellRanger (10× Genomics) with reads mapped to the human genome assembly (GRCh38).

We performed a secondary analysis on filtered matrix files using Seurat (v 3.2.0) (Stuart et al, 2019). Data were filtered for dead/poor quality cells based on the low number of genes detected (< 500) or potential doublets (> 6,000). Cells with a high percentage of mitochondrial genes were also removed at a threshold of < 35%. These filtered datasets were then scaled and normalized using SCTransform individually before being integrated based on 3,000 integration features. Clusters were identified in the integrated dataset using the FindCluster function at resolution 0.3, after PCA analysis, RunUMAP, and FindNeighbours at 1:30 dimensions. Cell type annotation was achieved through comparisons of various iterations from celldex (Aran et al, 2019) (R package) and known canonical markers of PBMC cell types. We finally settled on DICE annotation (Schmiedel et al, 2018). Differential expression analysis between NAFLD and healthy datasets (RNA) was performed using FindMarkers with MAST (Finak et al, 2015) (R package) for each individual cell type. Gene enrichment analysis was carried out for differential expression genes between PCV and normal in each cluster using clusterProfiler (R package, v 3.17.0.) (Yu et al, 2012).

Immune characterization by flow cytometry

Immune subsets were separated through MACS (Miltenyi Biotec). Phenotyping of both CD8+/− cell fractions was performed using flow cytometry antibodies directed against CD8‐PE (clone SK1, Biolegend 344706), CD4‐APC (clone RPA‐T4, Biolegend 300537). Details on antibodies are included in Appendix Table S5. Following 30 min of antibody incubation at room temperature, stained cells were washed and suspended in PBS containing 2% heat‐inactivated FBS. Flow cytometry was performed on BD LSRFortessaTM X‐20 Flow Cytometer (BD Biosciences). The gating was set by forward and side scatter, and 10,000 events were acquired. Finally, data were collected and analyzed using FlowJo v10.7.1 software (Becton Dickinson).

Endothelial‐leukocyte adhesion assay

The adhesion assay was as per described (Wilhelmsen et al, 2013). Briefly, THP‐1 and Jurkat cells were stained with a Deep Red CellTracker fluorescent dye (Invitrogen). 2.0 × 105 cells were then added to every cm2 of BOEC monolayers for 1 h. Subsequently, nonadherent cells were removed, and the number of adherent cells was determined through fluorescence intensity readouts from a Synergy H1 microplate reader (BioTek).

Endothelial‐leukocyte chemotaxis assay

CD8+/− immune cell fractions were isolated by magnetic‐activated cell sorting (130‐045‐201, CD8 MicroBeads, Miltenyi Biotec). Immune cells were allowed to rest for 24 h, then labeled with a Deep Red CellTracker (Invitrogen) prior to chemotaxis assay. The assays were performed on 96‐well Transwell plates with a membrane pore size of 5 µm, as per the manufacturer’s protocols (CBA‐105, Cell BioLabs). Briefly, BOEC monolayers were cultivated in EGM‐2 with 2% FBS for 48 h in the lower chamber before the CD8+/− cell fractions from the top chamber were allowed to migrate for 4 h. Migrated cells were then dislodged and collected from the underside of transwell membranes. The numbers of migrated cells were determined using ZEN 2 software (Zeiss).

Transendothelial migration assay

Blood outgrowth endothelial cells derived from three NAFLD subjects were pooled and seeded to confluence on the membranes with pores of size 5 µm (Corning) of the upper chambers of transwells for 24 h. The cells were then treated with or without LPS and FFAs for another 48 h. In parallel, PBMCs were isolated from blood samples of NAFLD subjects and allowed to rest in RPMI 1640 supplemented with 16% FBS for 24 h. On the day of assay, the PBMCs were resuspended in assay buffer (EGM‐2 + RPMI 1640 1:1 with 2% FBS and 0.5% BSA) and added to the BOEC monolayers in the upper chambers, and the PBMCs were allowed to transmigrate for 16 h. BOEC barrier integrity was characterized by transendothelial electrical resistance (TEER) by averaging probe measurements of three distributed spots over every transwell.

Inhibition of CXCL12‐CXCR4 axis in chemotaxis by chemokine receptor antagonist and chemokine neutralizing antibody

Plerixafor (AMD3100), a CXCR4 antagonist, was reconstituted as per the manufacturer’s instructions (HY‐10046, MedChemExpress). Immune cell fractions were then preincubated with 10 µM of AMD3100 for 30 min (Biasci et al, 2020). For CXCL12 inhibition, the BOEC monolayers in the lower chambers, and immune cells in the upper chambers were per‐incubated with an CXCL12‐neutralizing antibody (25 μg/ml, MAB310, R&D Systems) for 90 min (Zheng et al, 2020). Chemotaxis assays were then performed as per described above.

Profiling of circulating endothelial cells in NAFLD patients and healthy subjects

100 µl of 1 million PBMCs (whole PBMC fraction from processed blood) per subject were stained in the dark for 10 min at room temperature, followed by 20 min at 4°C on an analog tube rotator with antibodies (Appendix Table S5: key resource table). After incubation, cells were washed and resuspended in 200 µl of PBS with 1% BSA for flow cytometry analysis. CECs were detected through the combined immunophenotypic profile of CD45−/CD31+/CD133−/DNA+. Details on antibodies are included in Appendix Table S5. The number of CECs was expressed as cells per million of PBMCs. Flow Cytometry was performed using BD LSRFortessa X‐20 (BD Biosciences) and FACSDiva software (BD Biosciences), and data were analyzed using FlowJo v10.7.1 software (Becton Dickinson). Each analysis included at least 30,000 cells per individual.

Statistical analysis

Statistical significance of differences between the cohorts was analyzed using GraphPad Prism version 9. Data normality was determined either by Shapiro‐Wilk or Kolmogorov‐Smirnov test. Datasets with normal distributions were analyzed with the unpaired Student’s two‐tailed t‐tests to compare two conditions; one‐ or two‐way analysis of variance (ANOVA) followed by post hoc Tukey for datasets with more than two conditions. Nonparametric Mann–Whitney t‐test or Kruskal‐Wallis test were used for non‐normally distributed data. The application of these statistical methods to specific experiments is noted in the figure legends. A P‐value of < 0.05 was considered significant. Results are depicted as either mean ± SD or box plots indicating median (middle line), 25th, 75th percentile (box), and the lowest (respectively, highest) data points (whiskers).

Availability of materials

Further information and requests for resources and reagents should be directed to and will be fulfilled by the corresponding author. Some materials used in this study are commercially procured. There are restrictions to the availability of blood outgrowth endothelial cell lines derived from human patients and normal donors due to ethics considerations for use of these materials within the current scope of the study. Requests can be made to the corresponding author as we will explore the use of materials subject to new ethics approval and research collaboration agreement (including material transfer).

Author contributions

Chun‐Yi Ng: Data curation; Formal analysis; Validation; Investigation; Visualization; Methodology; Writing — original draft; Writing — review and editing. Khang Leng Lee: Data curation; Formal analysis; Validation; Investigation; Visualization; Methodology; Writing — original draft; Writing — review and editing. Mark Dhinesh Muthiah: Resources; Data curation; Formal analysis; Writing — review and editing. Kan Xing Wu: Data curation; Formal analysis; Investigation; Visualization; Methodology; Writing — review and editing. Florence Wen Jing Chioh: Data curation; Formal analysis; Visualization; Methodology; Writing — review and editing. Konstanze Tan: Data curation; Formal analysis; Visualization; Methodology; Writing — review and editing. Gwyneth Shook Ting Soon: Data curation; Formal analysis; Investigation; Writing — review and editing. Asim Shabbir: Resources; Writing — review and editing. Wai Mun Loo: Resources; Writing — review and editing. Zun Siong Low: Methodology; Writing — review and editing. Qingfeng Chen: Resources; Supervision; Methodology; Writing — review and editing. Nguan Soon Tan: Resources; Supervision; Writing — review and editing. Huck Hui Ng: Resources; Supervision; Funding acquisition; Project administration; Writing — review and editing. Yock Young Dan: Resources; Supervision; Funding acquisition; Project administration; Writing — review and editing. Christine Cheung: Conceptualization; Resources; Formal analysis; Supervision; Funding acquisition; Investigation; Visualization; Methodology; Writing — original draft; Project administration; Writing — review and editing.

In addition to the CRediT author contributions listed above, the contributions in detail are:

CC involved in conceptualization. C‐YN, KLL, MDM, KXW, FWJC, KT, and GSTS involved in data curation. C‐YN, KLL, MDM, KXW, FWJC, and CC involved in formal analysis. C‐YN, KLL, KXW, FWJC, KT, and GSTS involved in investigation. C‐YN, KLL, KXW, FWJC, KT, and ZSL involved in methodology. MDM, AS, WML, QC, and NST involved in resources. C‐YN and KLL involved in validation. C‐YN, KLL, KXW, FWJC, and KT involved in visualization. CC, YYD, and HHN involved in funding acquisition. CC, YYD, and HHN involved in project administration. CC, YYD, HHN, QC, and NST involved in supervision. CC, C‐YN, and KLL involved in writing—original draft. All authors involved in writing—review and editing.

Disclosure and competing interests statement

The authors declare that they have no conflict of interest.

Appendix

Click here for additional data file.

Expanded View Figures PDF

Click here for additional data file.