Inflammatory Alterations to Renal Lymphatic Endothelial Cell Gene Expression in Mouse Models of Hypertension

Hypertension (HTN) affects nearly half of all adults and leads to heart disease, stroke, and other severe cardiovascular conditions [1, 2]. The kidneys are actively involved in blood pressure (BP) regulation and are especially susceptible to HTN-related organ damage since they possess sensitive microvasculature that can be easily damaged by elevated BP, leading to immune cell activation and accumulation and chronic inflammation. As this process continues, pro-inflammatory and cytotoxic immune cells begin producing additional cytokines and attacking healthy cells [3, 4]. This process creates a loop of damage and inflammation that propagates the original issue of heightened BP. A great deal of work has been done to determine the contributions of immune cells in the kidneys to the development and maintenance of HTN. Notable experiments involve the depletion or adoptive transfer of immune cell subsets to and from hypertensive and normotensive animals, with several immune cell types taken from hypertensive animals being able to promote a pro-hypertensive state in normotensive animals [5, 6].

A natural control mechanism for this self-damaging feedback loop exists in the form of lymphatic endothelial cells (LECs), which are capable of trafficking or inactivating immune cells to reduce their inflammatory and cytotoxic potential in the surrounding tissue. Our lab has previously demonstrated that renal lymphatic vessels expand and grow as a result of HTN and hypertensive stimuli and that augmenting lymphangiogenesis in the kidney can reduce BP or prevent the development of HTN [7‒10]. However, little is known about the specific transcriptional changes that renal LECs undergo in response to HTN or how these changes enable them to compensate for elevated BP.

Single-cell RNA sequencing (scRNAseq) is an increasingly popular technique for providing an immense amount of transcriptional data in a variety of settings with incredible resolution [11]. This degree of resolution allows for a more refined understanding of transcriptomic changes and how specific cell populations can alter more granularly, rather than all at once. Additionally, new methods for analyzing these patterns and alterations allow for highly accurate predictions of cellular functionality that go beyond basic differential expression [12, 13].

In order to better evaluate the activity and response of renal LECs in HTN, we utilized scRNAseq on renal cell populations taken from hypertensive mice using two different mouse models of HTN. To combat the relative rarity of LECs in the kidney, we used a combination of the pan-endothelial marker CD31 and the lymphatic marker podoplanin, which are common markers to locate and specify renal LECs and lymphatic vessels [14, 15]. Single-cell transcriptomics allow for an exceptional degree of resolution in determining how individual cell types respond on a functional level to a condition and we sought to use this approach to better dissect how renal LECs respond to the murine angiotensin II-induced (A2HTN) and salt-sensitive (SSHTN) models of HTN.

Male C57/BL6 mice were purchased from Jackson Laboratories and were between 10 and 14 weeks of age at the start of each model. We utilized the A2HTN and SSHTN models as described previously (n = 6 for each group) [7, 9]. Briefly, for the A2HTN model, mice were implanted with subcutaneous osmotic pumps under anesthesia. The pumps delivered either 1,000 ng/kg/min of angiotensin II in saline or pure saline for 3 weeks prior to euthanasia by isoflurane followed by organ retrieval. For the SSHTN model, mice were primed with 0.5 g/L of LNAME in their drinking water for 2 weeks followed by a 2-week washout period. Afterward, mice received either 4% salt chow or normal chow diet for 3 weeks prior to euthanasia and organ retrieval. Animal use protocols were approved by the IACUC of Texas A&M University.

Systolic BPs were measured by tail cuff using the IITC Life Science BP acquisition system as described previously [16]. The chamber and restraints were warmed to 34°C and mice were allowed to acclimate for 5 min prior to measurement. BPs were taken prior to surgery and on days 7, 14, and 21 for the A2HTN groups or at the end of the washout period and on days 7, 14, and 21 during the 4% salt diet for the SSHTN groups.

At the end of each model, the animals were euthanized and the kidneys from all animals in a given group were pooled and finely minced. Enzymatic digestion was performed as described previously using a combination of Collagenase D (2.5 mg/mL, Roche) and Dispase II (1 mg/mL, Sigma-Aldrich) and was accompanied by mechanical dissociation through light shaking at 10-min intervals [7, 17]. The cell suspension was filtered using 70 μm and 40 μm filters and washed twice to remove dead cells. Cells were then incubated with biotinylated rat anti-mouse CD31 antibodies (BD Biosciences, catalog number: 553371) and separated using the Invitrogen CELLection Biotin Binder Kit (Thermo Fisher Scientific, catalog number: 11533D). This process was repeated using biotinylated goat anti-mouse podoplanin antibodies (R&D Systems, catalog number: BAF3244) and the Invitrogen CELLection Biotin Binder Kit, then filtered and washed again. Cells were counted and assessed for viability using the Vi-CELL XR Cell Counter (Beckman Coulter Life Sciences).

Libraries were prepared with the 10X Genomics Chromium Next GEM Single Cell 3′ Reagent Kit V3.1 using its standard protocol and cDNA from each sample was verified with the Agilent High Sensitivity D5000 ScreenTape assay, quantified, normalized, and then sequenced on an Illumina NovaSeq 6000 with an S4 flow cell by the Texas A&M Institute for Genome Sciences and Society. The resulting FASTQ files were aligned with Cell Ranger 7.0.0 via 10X Genomics Cloud Analysis using the 10 mm (2020 revision) mouse genome assembly.

The files for each group were loaded into Seurat version 4.1 [18]. For quality control, genes and cells were removed from each sample using parameters derived from previous work in scRNAseq analysis of renal endothelium [19]. Cells with <300 or >3,000 genes and cells with a mitochondrial gene ratio above 10% were removed and genes expressed by <50 cells and mitochondrial genes were removed. Additionally, DoubletFinder 2.0.3 was run to exclude likely doublets in the sample [20]. Cells were clustered using 20 principle components for PCA and each cluster’s cell type identity was determined based on a list of genes derived from bulk RNAseq and scRNAseq data of renal epithelium, blood and lymphatic endothelium, and immune cells [21‒23]. Based on the expression and/or absence of known markers and their 100 most highly expressed genes, we identified groups of LECs, myeloid immune cells, and an unexpected group of multipotent cells, then labeled the clusters accordingly.

Due to the relatively small number of LECs in the control samples and the lack of significant differences in the BP and gene expression between the two models’ controls, the two control samples were merged and reclustered into a single control group to improve clustering resolution for the purpose of comparative analyses. Differential expression between hypertensive (HTN) and normotensive control (Ctrl) samples was performed in Seurat with default parameters. Gene set enrichment analysis (GSEA) was performed using Fast Gene Set Enrichment Analysis (FGSEA 1.20) and the mouse curated and ontology gene sets from the Molecular Signatures Database via msigdbr 7.5.1 (https://www.gsea-msigdb.org/gsea/msigdb/) [24].

The Python implementation of SCENIC (pySCENIC 0.11.2 in Python 3.8) and mouse transcription factor motifs (v9) were used to find differentially active positive regulons [25]. Regulon prediction was performed using default parameters and the resulting matrices were used to determine regulon specificity scores between HTN and Ctrl groups. Regulons are referred to in the text, figures, and tables by the identified genes that regulate the full group of enriched genes.

Weekly tail cuff BP readings after the start of each model (days 7, 14, and 21) were elevated in HTN groups over their respective controls (Table 1; p < 0.001), similar to previously reported values [7, 9, 10]. The A2HTN control group (A2 Ctrl) and SSHTN control group (SS Ctrl) showed no significant differences in BP, which aligns with previous data showing minimal differences between these control groups.

Samples were grouped into merged HTN and merged Ctrl aggregate groups to determine model-neutral effects of HTN. After FASTQ alignment, the samples were loaded into Seurat for quality control for cell health and doublets, then the original and merged groups were clustered and cluster-specific gene markers were analyzed [18]. A total of ∼15,000 cells were present across all samples and several distinct and consistent clusters were identified in each sample and indicated the presence of multiple cell types in addition to multiple cellular states per cell type.

Considering the dual enrichment for CD31+ and podoplanin+ and the variety of known cell types present in the kidney, as well as the markers they express as mRNA and cell surface markers, we determined that there were renal LECs, myeloid immune cells (MICs), and a third, unusual cell type which we dubbed support cells (SCs) present in these samples. Clusters in each sample were labeled as one of these cell types using known cell type-specific markers from bulk RNAseq and scRNAseq (Fig. 1 and online suppl. Fig. S1, S2; for all online suppl. material, see https://doi.org/10.1159/000539721; Tables 2, 3 and online suppl. Tables S1–S7) [19, 21‒23, 26‒41]. We noted a relative scarcity of LECs in the SS Ctrl sample, which had led to issues with clustering resolution, and chose to utilize the merged Ctrl group for all further comparisons involving a control sample to avoid this issue and the potential bias it could introduce. The top 200 most variable genes per cell type were determined across all samples to help verify the assigned labels (Fig. 2; online suppl. Table S8).

To focus on the effects of HTN on LECs, we isolated and reclustered the LECs present in each group then reanalyzed the cluster-specific markers (Fig. 3 and online suppl. Fig. S3–S4; Tables 4, 5 and online suppl. Table S9–S13) [42‒49]. Broadly speaking, there was a marked increase in the number of LECs present in the HTN samples over the Ctrl samples as well as an increase in the heterogeneity of cellular states, which aligns with known pathological alterations in LECs [8, 50‒53].

Differentially expressed genes (DEGs) were identified with Seurat using the merged control LEC group and an adjusted p value of <0.01 for all comparisons. When comparing all HTN and all Ctrl LECs, there were a total of 597 DEGs (486 ↑, 112 ↓; online suppl. Table S14). Among the top 25 most significant DEGs are genes related to growth, inflammation, and cytoskeletal remodeling (Table 6) [42, 43, 54‒69].

In the top 25 DEGs, CD31 (Pecam1) and neuropilin-1 (Nrp1) were upregulated by 69% and 67% respectively and both have functions related to both endothelial growth and immunoregulation [70, 71]. In addition, lymphocyte antigen 6 family member E (Ly6e), a gene known to enhance immune checkpoint activity as well the growth-related effects of hypoxia-inducible factor-1 alpha subunit (HIF-1α; Hif1a), was also significantly upregulated with an 81% increase in expression and was expressed in 95% of LECs in HTN (up from 65% in Ctrl) [72, 73]. Other immune checkpoint and immunomodulatory genes were also upregulated in HTN, including butyrophilin-like 9 (Btnl9) (86% increase), TRAIL (Tnfsf10) (42% increase), galectin 9 (Lgals9) (31% increase), and CD47 (26% increase) [74, 75]. Similar to a previous report involving the effects of Leishmania major infection, we also noted an increase in H2-K1, H2-D1, and beta-2-microglobulin (B2m), which are all a part of MHC class I antigen presentation in LECs, as well as increases in guanylate-binding protein (GBP) genes, which is induced by interferon signaling and can lead to inflammasome activation [76‒78].

Comparisons for each individual model were also conducted, using A2HTN or SSHTN LECs instead of the merged HTN LEC group. A2HTN LECs had a total of 1,070 DEGs (948 ↑, 122 ↓; online suppl. Table S15) whereas SSHTN LECs only had 315 DEGs (244 ↑, 71 ↓; online suppl. Table S16). The DEGs for each model were largely similar to the merged HTN comparison, though Ly6e’s upregulation was even more exaggerated in the A2HTN LECs and was accompanied by upregulation of additional immune checkpoint genes, including PD-L1 (Cd274) (24% increase), CD200 (30% increase), nectin cell adhesion molecule 2 (Nectin2) (30% increase), and B2m (87% increase) [75, 79‒81]. Adding to the upregulation of MHC I-related genes seen in the merged HTN comparison, there was upregulation of MHC class II genes (notably CD74, H2-Ab1, and H2-Aa) in A2HTN LECs. In SSHTN, heat shock protein 90 alpha family class B member 1 (Hsp90ab1) increased by 93% (vs. a 22% increase in A2HTN) and heat shock protein family B (small) member 1 (Hspb1) increased by 114% (vs. 32% in A2HTN), possibly indicating a preference for heat shock protein signaling for growth in the presence of both excessive BP and sodium [54, 57, 60, 82]. When comparing LECs from the two models directly, additional growth and development genes are upregulated in A2HTN over SSHTN, including CD151, which has been associated with lymphangiogenesis and invasion in cancer, and B2m and Ly6e are increased by 48% and 32% respectively (online suppl. Table S17) [83].

GSEA was performed using FGSEA with the identified DEGs and the curated gene sets (M2) and ontology gene sets (M5) from the mouse collections of the Molecular Signatures Database. Pathways were filtered to only contain entries with adjusted p values <0.01 and are shown sorted by normalized enrichment score (NES) in figures. Comparing all hypertensive LECs to normotensive LECs, we identified 56 significantly enriched pathways (18 ↑, 38 ↓) for the M2 gene sets and 218 (62 ↑, 156 ↓) for the M5 gene sets (Fig. 4; online suppl. Tables S18, S19). A2HTN LECs compared to normotensive LECs yielded 75 (35 ↑, 40 ↓) and 181 (51 ↑, 130 ↓) enriched pathways for the M2 and M5 sets respectively, while SSHTN LECs had 58 (14 ↑, 44 ↓) and 227 (58 ↑, 169 ↓) pathways enriched over Ctrl LECs (online suppl. Fig. S5; online suppl. Tables S20–S23). 19 (17 ↑ in A2HTN, 2 ↑ in SSHTN) and 12 (10 ↑ in A2HTN, 2 ↑ in SSHTN) pathways for M2 and M5 sets were differentially enriched between the two HTN models (online suppl. Tables S24–S25). Overall, the enriched pathways indicate a strong inflammatory and immune-type response, the induction of robust angiogenic/lymphangiogenic growth mechanisms, and an overall shift in metabolic processes in HTN. A2HTN LECs again show an exaggerated inflammatory response in addition to an increase in immune markers whereas SSHTN LECs seem to be more characterized by the downregulated pathways, as all pathways common to the two models were more downregulated in SSHTN LECs than in their A2HTN counterparts.

Utilizing the resolution provided by scRNAseq, new approaches have been developed to identify regulons, which consist of a group of genes which are expressed as a unit and under the control of a single master gene (usually a transcription factor). By evaluating co-expression of experimentally-validated master genes and groups of their subordinate genes in different conditions, differential regulon enrichment can be identified and help provide more detailed context to known alterations in transcription factor activity and how they relate to the broader pathophysiology of a disease, including how pathologically altered genes influence one another. To determine if any gene expression changes in LECs could be tied to transcription factor activity related to HTN, regulons and regulon specificity scores (RSSs), which measure how specific a regulon’s activity is to a given group, cluster, or condition, were identified using pySCENIC for each group and each comparison respectively [25, 84].

RSS values were determined for merged hypertensive, A2HTN, and SSHTN LECs versus normotensive LECs, as well as comparing LECs from each model (online suppl. Tables S26–S29). A total of 605 regulons were identified for the HTN versus Ctrl LEC comparison, with 331 having RSS values >0.5 for HTN LECs, indicating preference for activation in HTN, and 64 having RSS values >0.5 for Ctrl LECs, indicating preference for activation in controls. A2HTN versus Ctrl LECs showed 602 regulons (342 for A2HTN, 72 for Ctrl) and SSHTN versus Ctrl LECs gave 583 regulons (345 for SSHTN, 57 for Ctrl). Comparing the two HTN models yielded 613 regulons, with 192 showing preference for A2HTN and 263 with preference for SSHTN. The regulons with the highest RSSs for each hypertensive versus normotensive comparison were plotted by expression and RSS (Fig. 5 and online suppl. Fig. S6) and clustermaps were generated for all regulons in each group (online suppl. Fig. S7).

The top five regulons by RSS for HTN LECs were MAF BZIP transcription factor F (Maff), high mobility group nucleosomal binding domain 3 (Hmgn3), HIV type I enhancer-binding protein zinc finger 1 (Hivep1), Ras responsive element binding protein 1 (Rreb1), and Hif1a. Rreb1 and Hif1a are both major drivers for growth in LECs in inflammatory conditions, while Hivep1 negatively regulates NF-κB signaling in monocytes and macrophages [85‒90]. Additionally, Maff has anti-inflammatory functions and can alter the function of neighboring monocytes and Hmgn3 is upregulated by TGFB1 and the high-mobility group (HMG) family of proteins, which includes Hmgn3, has a strong association with immune checkpoint genes [91‒97].

The A2HTN and SSHTN models used here generated robust increases in systolic BP which align with previous studies from our lab and others and are known to be sufficient to cause serious renal vascular injury [9, 98, 99]. The cell surface marker enrichment for CD31+/podoplanin+ cells was successful in increasing the proportion of LECs present in the sample, but it also had the effect of pulling in MICs and SCs [14, 15, 100]. While unexpected, these cell types are both known to express these markers under different conditions and may represent distinct subsets of larger populations that weakly express similar markers to LECs or even participate in inflammatory processes akin to LECs. The SCs we identified share many characteristics with mesenchymal stem cells, which are known to have anti-inflammatory effects and can even reduce BP and improve kidney function, while MICs can have a wide variety of pro-inflammatory, anti-inflammatory, and pro-lymphangiogenic effects [6, 101‒104].

Our previous studies have shown an increase in LECs and LEC markers in the kidneys in HTN and our results here seem to reinforce this at the single-cell expression level. Not only do we see an increase in LECs and genes known to be involved in lymphangiogenesis in the HTN samples, but we see an increase in the diversity of transcription factors and genes expressed by LECs, though some LEC markers appear to differ at single-cell resolution and across different organs [22, 105, 106]. However, while there were multiple clear LEC clusters in the HTN samples, we had far fewer cells that could be grouped into obvious LEC clusters in control samples, which likely reflects a relative scarcity of LECs in kidneys under normal conditions. Additionally, most DEGs when comparing hypertensive and normotensive LECs were upregulated rather than downregulated (ranging from 77 to 89% of all DEGs being upregulated in HTN, depending on the comparison), possibly indicating the LECs exist in a fairly quiescent state under normal conditions, then react strongly in response to inflammation.

Regarding canonical LEC and lymphangiogenic markers in these samples, we observed that SRY-box transcription factor 18 (Sox18) is upregulated by 65% in merged HTN LECs and is up by 105% in A2HTN LECs, but does not have any sort of regulon specificity between the groups. On the other hand, COUP-TF II (Nr2f2) is not differentially expressed, but has specificity scores for HTN groups ranging from 0.721 to 0.767. Similarly, the master lymphatic transcription factor prospero homeobox 1 (Prox1) is not strongly or differentially expressed in these samples, yet it is among the top variable genes for LECs in this dataset and has a slightly higher RSS for A2HTN compared to controls (0.250 vs. 0.167), and VEGFC also ranks in the top 100 variable genes for both LECs and MICs.

Conversely, while many of the DEGs, pathways, and regulons presented here are related to inflammation, growth, and expansion in endothelial cells, they are not a complete match for the traditional suite of genes associated with lymphangiogenesis. Instead, many of them are more comparable to studies involving various forms of inflammatory processes in LECs in other conditions or states and may indicate that these new lymphatics are being generating by and expanding through a continuous flood of HIF-1α, VEGF, and TLR signaling, rather than the more coordinated signaling that triggers and sustains lymphangiogenesis during development [73, 107‒110]. Notably, HIF-1α and its downstream effectors can contribute to lymphangiogenesis and lymphatic vessel malformation when upregulated [88]. PI3K (Pik3ca), which is one of HIF-1α’s effectors, has a well-known mutation in cancers that exaggerates lymphangiogenesis and lymphatic malformations, with recent scRNAseq analysis of this mutation’s activity revealing a number of similarities to the genes and pathways we have presented here [111].

In addition to growth-related effects, HIF-1α can also directly stimulate the expression of immune checkpoint ligands, including PD-L1, which can regulate and reduce the number of cytotoxic T cells in an area [85]. Immune checkpoint inhibitors as part of cancer treatment have been observed to have negative effects on systemic and pulmonary HTN, though this is somewhat disputed [112‒114]. However, the relationship between HTN and activation of the immune system is well-established and is a topic of ongoing research, particularly regarding how the immune system can directly impact BP [115, 116]. CD8+ T cells are involved in the onset and worsening of HTN and other cardiovascular diseases through their pro-inflammatory and cytotoxic functions, so reducing their numbers and activity in the renal interstitium through lymphatic immune checkpoint signaling would help mitigate further inflammation and damage, potentially even lowering BP as a result [3, 4, 117]. Through the upregulation of genes involved in immune checkpoints and immunoregulation, including Ly6e, Btnl9, Lgals9, B2m, Tnfsf10, and Cd274, renal LECs (particularly in the A2HTN model) appear to be increasing their immunosuppressive potential, which is likely a natural compensatory response to limit the excessive and maladaptive inflammation that is present in the kidneys in HTN.

To better understand how renal LECs are altered by and influence the progression of HTN, we performed scRNAseq on CD31+/podoplanin+ renal cells taken from A2HTN and SSHTN mice and their respective controls. We identified LEC, SC, and MIC populations and analyzed the changes in gene expression in LECs due to HTN through differential gene expression, GSEA, and regulon analysis. There was a marked increase in the number and activity of LECs in the HTN groups, including the upregulation of inflammatory growth pathways that are similar to what is seen in other inflammatory conditions. Additionally, multiple genes related to immune checkpoint and immunoregulatory signaling were upregulated in LECs in HTN, providing a new potential role for renal LECs in counteracting the development and maintenance of HTN and HTN-related kidney damage.

In the future, we plan to evaluate how the inflammatory growth seen here impacts the function of LECs and their role in immune cell trafficking. As this type of growth has been associated with lymphatic hyperplasia and malformations, there may be a substantial impact on clearance and trafficking, leading to an accumulation of pro-inflammatory immune cells in the kidneys. Also, with the increase in immune checkpoint genes, we will investigate how lymphatic immune checkpoint and immunoregulatory signaling regulates the immune cell populations responsible for the onset and continuation of HTN. Better understanding of the pathological changes in renal LECs and immune cells that influence the development and progression of HTN will help provide novel approaches for treating severe HTN.

This study protocol was reviewed and approved by the IACUC of Texas A&M University, Approval No. 2022-0083.

The authors have no conflicts of interest to declare.

This work was funded by NIH RO1 (DK120493) to B.M. Mitchell and Texas A&M University T3 to J.M. Rutkowski.

Justin G. McDermott: conceptualization, data curation, formal analysis, investigation, methodology, project administration, resources, software, writing – original draft, and writing – review and editing. Bethany L. Goodlett and Shobana Navaneethabalakrishnan: data curation and investigation. Heidi A. Creed: conceptualization, data curation, and investigation. Joseph M. Rutkowski: conceptualization, funding acquisition, and project administration. Brett M. Mitchell: conceptualization, funding acquisition, project administration, resources, and writing – review and editing.

The data that support the findings of this study are publicly available in NCBI GEO (GSE 236410) at https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE236410. Any additional inquiries can be addressed by the corresponding author.