Aortic heterogeneity across segments and under high fat/salt/glucose conditions at the single-cell level

Abstract The aorta, with ascending, arch, thoracic and abdominal segments, responds to the heartbeat, senses metabolites and distributes blood to all parts of the body. However, the heterogeneity across aortic segments and how metabolic pathologies change it are not known. Here, a total of 216 612 individual cells from the ascending aorta, aortic arch, and thoracic and abdominal segments of mouse aortas under normal conditions or with high blood glucose levels, high dietary salt, or high fat intake were profiled using single-cell RNA sequencing. We generated a compendium of 10 distinct cell types, mainly endothelial (EC), smooth muscle (SMC), stromal and immune cells. The distributions of the different cells and their intercommunication were influenced by the hemodynamic microenvironment across anatomical segments, and the spatial heterogeneity of ECs and SMCs may contribute to differential vascular dilation and constriction that were measured by wire myography. Importantly, the composition of aortic cells, their gene expression profiles and their regulatory intercellular networks broadly changed in response to high fat/salt/glucose conditions. Notably, the abdominal aorta showed the most dramatic changes in cellular composition, particularly involving ECs, fibroblasts and myeloid cells with cardiovascular risk factor-related regulons and gene expression networks. Our study elucidates the nature and range of aortic cell diversity, with implications for the treatment of metabolic pathologies.

4℃ . After reverse transcription, GEMs were broken and the single-strand cDNA was cleaned up with DynaBeads MyOne Silane Beads (Thermo Fisher Scientific, USA). cDNA was amplified on an S1000 TM thermal cycler (Thermo Fisher Scientific) running the following program: 98℃ for 3 min, 98℃ for 15 s, 67℃ for 20 s, 72℃ for 60 s, 12 cycles; 72℃ for 1 min, and held at 4℃ . After cDNA amplification, cDNA was cleaned with SPRIselect Reagent (Beckman Coulter) and the cDNA sample was loaded onto an Agilent Bioanalyzer High Sensitivity chip for qualitative analysis. Indexed sequencing libraries were constructed using the reagents in the Chromium TM Single Cell 3' Reagent Kits (10x Genomics) in the following steps: (1) fragmentation, end repair, and A-tailing; (2) end repair and A-tailed double-sided size selection with SPRIselect reagent; (3) adaptor ligation with an Adaptor Ligation Mix and then post-ligation cleanup with SPRIselect reagent; and (4) sample index PCR and cleanup. The Agilent Bioanalyzer High Sensitivity chip was used for qualitative analysis of the library.
Sequencing was performed on an Illumina HiSeq X Ten platform. The accession number for the RNA-seq data reported in this paper is NCBI SRA: PRJNA489757.
Sequencing data pre-processing CellRanger software (10x Genomics) was then used to analyze the sequencing data and produce gene expression information for each cell 84 . The software converted Illumina basecall files to fastq format. The clear raw sequence data (fastq) was produced using the FASTX-Toolkit. Then the reads were aligned to the mouse mm10 transcriptome. The total UMI counts for each cell-associated barcode were normalized to the grand median UMI count per cell by a scaling factor (median reads per barcode / reads per barcode), and then principal component analysis (PCA) was used to reduce the dimensionality of the raw data. The output dataset was loaded into Cell Ranger (10x Genomics) to generate graph-based and K-mean clusters in the PCA space. In this study, a graph-based algorithm was used in dimensionality-reduction to generate two-dimensional t-distributed stochastic neighbor embedding (t-SNE) plots. Single-cell expression data, including log 2 fold-change (FC) values (negative binomial test) and mean UMIs (normalized by the size factor: total UMI count per cell / median UMI count per cell) were generated with CellRanger. For quality control, cells with >10% reads mapping to mitochondria, or with <700 or >15,000 UMI counts indicated low quality, and were removed from the analysis.
Wire myography Vascular reactivity of freshly isolated aortic arteries was studied using a myograph (Danish Myo Technology, Aarhus, Denmark) as described 86 . After mice were sacrificed, the aorta was rapidly removed and placed in oxygenated ice-cold Krebs solution that contained (mmol/L): 119 NaCl, 4.7 KCl, 2.5 CaCl 2 , 1 MgCl 2 , 25 NaHCO 3 , 1.2 KH 2 PO 4 , and 11 D-glucose. Connective tissue and blood were removed and four segments (2 mm in length) were mounted in organ baths individually gassed with 95% O 2 and 5% CO 2 and maintained at 37°C. Changes in isometric tone of the rings were recorded. The rings were stretched to an optimal baseline tension of 3 mN. After 60 min of equilibration at baseline tension, 60 mM KCl (the NaCl in Krebs solution was substituted with an equimolar amount of KCl) was applied to confirm the contractility of the ring, washed in Krebs solution, and finally allowed to equilibrate for 30 min. Phenylephrine (Phe) (1 nmol/L to 1 μmol/L) was added cumulatively, or Phe (1 μmol/L) was used to produce a steady contraction and then acetylcholine (ACh) (1 nmol/L-10 μmol/L) was added cumulatively. The dose-dependent contractile response to Phe is expressed as a percentage of 60 mmol/L KCl-induced contraction, while relaxation responses to ACh are expressed as a percentage of the constriction induced by Phe (1 μmol/L).
Evans blue permeability assay Mice were anesthetized with 4% chloral hydrate and fixed on a plate. Evans blue dye (5 mg/kg 1% Sigma) was injected into the inferior vena cava to stain the aorta through the systemic circulation. After 30 min of incubation, the intact aorta was removed and photographed under a light microscope. To quantify the aortic permeability, free Evans blue was removed by PBS perfusion. Then the aorta was weighed and cut into pieces. The Evans blue in tissue was extracted with formamide overnight at 55°C and the concentration of Evans blue was measured on a spectrophotometer at 620 nm. The Evans blue content per gram of tissue was calculated.

STATISTICAL ANALYSIS
Cell clustering Principal component analysis was used to decompose the expression matrix of each cell into linear components by CellRanger software implemented with the IRLBA R package. Then the cells were clustered with a graph-based clustering algorithm and visualized as t-SNE plots using CellRanger software. The cell types were determined using significantly up-regulated marker genes (log 2 fold-change >1.5, p <0.05 versus the other clusters by the negative binomial test) in a certain cluster in the t-SNE plot, which was visualized with either Loupe Cell Browser (10x Genomics) or in the cell expression matrix generated by CellRanger.

GO analysis
To identify the potential functional pathway of a set of genes, GO enrichment analyses were performed using the DAVID 6.8 functional annotation tool 88 . Briefly, cell types and subpopulations were determined using the graph-based t-SNE plots. A gene expressed with a log 2 FC value >1.5 (negative binomial test) and a p value <0.05 between the indicated cell types/subpopulations and all remaining groups was determined to be up-regulated between the cell types/subpopulations. Then the list of significantly changed genes served as the input of GO in DAVID. A significantly enriched functional pathway was determined by the enrichment score and p value generated by DAVID.
Gene regulatory network analysis Regulon matrices were generated by SCENIC and used to map gene-regulatory networks 89 . The scRNA-Seq data from the abdominal aorta in healthy, high-salt intake, high-fat intake, and high plasma glucose mice were combined, and then ECs, myeloid cells, and fibroblasts were extracted from the combined dataset with known marker genes -PECAM1 for ECs, CD163 for myeloid cells, and PDGFRA for fibroblastsusing CellRanger software. Gene-expression matrices in which rows corresponded to log-transformed gene values and columns corresponded to cells were produced for each cell type. Then the matrices were loaded into GENIE3 to infer transcription factor co-expression modules 37 . Co-expression modules with <30 genes were filtered out. RcisTarget was then used for cis-regulatory motif analyses to build regulons. We found 329 regulons in ECs, 315 in fibroblasts, and 317 in myeloid cells. The activity of the regulons in cells was then calculated as the area under the recovery curve by AUCell across the genes ranked by their expression values, and active regulons were determined by the AUCell default threshold.
The binary regulon activity matrix was mapped to t-SNE plots with the Rtsne package.
Regulons IRF9, NR2C2, and ATF1 were chosen as they predominated in at-risk ECs, fibroblasts, and myeloid cells but were inactive in healthy cells.
To determine the potential for fibroblast-EC communication in the aorta, we used mouse orthologs of human ligand-receptor pairs designed by Ramilowski et al. from the BioMart database 57 . We defined one fibroblast or endothelial subpopulation as expressing a certain ligand or receptor if the gene was significantly higher in the subpopulation than in the other subpopulations with a log 2 FC value >1.5 (negative binomial test) and p <0.05. When fibroblasts and ECs were in communication, the ligand was expressed in the fibroblast subpopulation and the receptor was expressed in the EC subpopulation 80 . Then a ligand-receptor matrix of the cell identities, log 2 FCs, and p values of genes was created, and the list of receptors in the matrix was used for GO analysis. A heatmap (representative in are not included. Data are shown as the mean ± s.e.m. c. Contribution of batches to each cluster. Dot-plot and t-SNE showing the batch composition of each detected cluster crossing four aortic segments from controls or at-risk aortas. All replicates contributed to all clusters, and no major batch effect was found. A total of 32 data from all aortic segments were used. Data are shown as the mean ± s.e.m. d,