Thirteen adult male mice (C57BL/6J; 14–20 weeks old, Jackson Laboratory, Bar Harbor, ME, United States) were housed in an antigen-free and virus-free animal care facility under a 12-h light and dark cycle. Six mice were exposed to hypoxic conditions by placing them in a plexiglass chamber connected to OxyCycler Models A420OC and AT42CO gas controllers and analyzers (BioSpherix, NY, United States). This gas control delivery system regulated the flow of room air, N₂, and O₂ into the chamber to create and maintain prescribed fraction inspired O2 (FiO₂) levels. Persistent hypoxia defined by a FiO₂ of 10% was maintained for 3–6 weeks at 24 h/day. Seven age matched C57BL6J mice were kept in similar conditions but exposed to room air (∼FiO₂ 20%) as controls. Both groups of mice were fed a standard rodent chow and had free access to water. Animals were euthanized with overdose of urethane by intraperitoneal injection followed by exsanguination, followed by harvest of the hearts, lungs, and pulmonary arteries.

Specimens were excised from the main pulmonary artery to the first branch of the right pulmonary artery (RPA) and prepared as previously described (Ramachandra and Humphrey, 2019). After flushing of blood with a Hanks buffered physiologic solution (Hanks and Wallace, 1949), perivascular tissue and fat were gently removed, and the left pulmonary artery and small branch vessels were ligated with suture. The RPA was cannulated on custom glass micropipettes and secured beyond the main pulmonary artery with ligature on one end and the first branch of the pulmonary at the other end. The specimen was submerged in a 37°C bath of Krebs-Ringer solution perfused with 95% O₂/5% CO₂. The specimens were tested using a custom computer-controlled testing device (Gleason and Humphrey, 2004), with an active protocol followed by a passive protocol (n = 5 for hypoxia, n = 6 for normoxia). Active testing focused primarily on SMC contractile responses near in vivo conditions, that is, a mean distending pressure of 15 mmHg and the specimen-specific value of in vivo axial stretch as previously described (Murtada et al., 2016; Caulk et al., 2019; Ramachandra and Humphrey, 2019). Such protocols are much more relevant physiologically than the standard “ring-tests” performed with uniaxial myographs. First, however, the SMCs were conditioned at 10 mmHg and an axial stretch of 1.1 (relative to the unloaded length) by contracting the vessel with 100 mM KCl. Axial stretch was subsequently increased in steps of 0.01 and pressure in steps of 1 mmHg to avoid sudden large deformations, which could compromise the viability of the SMCs and ECs. Following conditioning, the vessels were contracted with 100 mM KCl and then 1 mM phenylephrine (PE), both at 15 mmHg and the in vivo axial stretch. This protocol consisted of 5 min of equilibration, 15 min of contraction, and 10 min relaxation following washout of the vaso-stimulant with a fresh Krebs-Ringer solution. EC testing immediately followed maximal SMC contraction with PE without vaso-stimulant washout. Specimens were exposed to 10 μM acetylcholine (ACh) followed by inhibition of eNOS release by 1 mM N(gamma)-nitro-L-arginine methyl ester (L-NAME). Arterial diameter is regulated by the combined activities of EC and SMC; therefore, inhibition of EC secretion of vasodilatory molecules such as nitric oxide enables determination of maximal SMC contractility (Murtada and Humphrey, 2018).

A Titanium-Sapphire Laser (Chameleon Vision II, Coherent) was used to image representative regions of pulmonary arteries at in vivo relevant loading conditions (in vivo stretches and pressures identical to those used during passive mechanical testing, described below). A LaVision Biotec TriMScope microscope was tuned at 840 nm and equipped with a water immersion 20× objective lens (NA. 0.95). The backward scattering second harmonic generation signal from fibrillar collagens was detected within the wavelength range 390–425 nm; the autofluorescent signal arising from elastin was detected at 500–550 nm, and the fluorescent signal of cell nuclei labeled with Syto red stain was detected above 550 nm. An in-plane field of view (axial-circumferential plane) of 500 μm × 500 μm and a volume of about 0.05 mm³ were used; this provides a much greater volume of tissue for imaging than via standard histology and hence averaging over significantly greater numbers of cells and extracellular matrix. The in-plane resolution was 0.48 μm/pixel and the out-of-plane (radial direction) step size was 1 μm/pixel. 3D images acquired concurrently for the three signals (collagen, elastin, and cell nuclei) were post-processed using MATLAB R2019b and ImageJ 1.53a. The first processing step relied on the near cylindrical shape of the samples to fit a circle to the two-dimensional mid-thickness profile of the arterial wall and transform each circumferential-radial slice of the 3D images from Cartesian to polar coordinates (angle and radius). This allowed a layer-specific microstructural analysis to focus on collagen fiber alignment and cell volume density analyses, as described in previously (Cavinato et al., 2020, 2021).

Following biaxial testing and multiphoton imaging, the specimens were fixed in 10% neutral buffered formalin and stored in 70% ethanol at 4°C for histology. After embedding in paraffin, they were sectioned (5 μm thickness) and radial-circumferential cross-sections (planes) were stained with Verhoeff Van Giesen (VVG), Masson’s Trichrome (MTC) or Movat pentachrome. Details of image quantification can be found elsewhere (Bersi et al., 2017). Briefly, each section was imaged with an Olympus BX/51 microscope using an Olympus DP70 digital camera under a 20× magnification objective. Complete cross-sections were obtained by stitching together sub-images with Image Composite Editor software (Microsoft Research). The stitched images were subsequently analyzed using custom MATLAB scripts. Briefly, following background subtraction and pixel-based thresholding, area fractions for elastin (from VVG) and cytoplasm (from MTC) were computed as the ratio of pixels corresponding to a stain divided by the total number of pixels in the image. Because MTC can overstain collagen, its area fraction was computed as 1 - area fraction of elastin plus cytoplasm, with GAG content assumed negligible (Ferruzzi et al., 2015, 2018) as confirmed with Movat Pentachrome staining (data not shown). Three sections were analyzed per vessel per stain.

We used a 2-D formulation to model the passive mechanical behavior since residual stresses tend to homogenize the stress field, thereby rendering mean values as good estimates of overall wall stress (Humphrey, 2013). Mean stresses along the circumferential (θ) and axial (z) directions are,

where P is the transmural pressure, r_i internal radius, r_o outer radius, and f_T the transducer measured axial force. Under the assumption of incompressibility, inner radius ri=ro2-V¯/πl, where l is the instantaneous length between the ligatures securing the vessel to the micropipettes and V¯ is the volume of the vessel in the unloaded state; V¯=πL(OD2-ID2)/4 where L is the unloaded length, OD the unloaded outer diameter, and ID = OD −2H, the unloaded inner diameter.

A pseudoelastic constitutive formulation, which has successfully described passive biaxial behaviors of pulmonary arteries (Ramachandra and Humphrey, 2019), modeled the passive mechanical behavior in terms of a stored energy function, W. Associated wall stress and material stiffness can be computed from first and second derivatives of W with respect to an appropriate deformation metric. Based on success of prior work, we let

where c, c1i and c2i are material parameters; i = 1, 2, 3, 4 represent four collagen-dominated families of fibers along axial, circumferential, and two symmetric diagonal directions, respectively. I_c is the first invariant of right Cauchy-Green tensor (=λθ2+λz2+1/λθ2λz2) and IVci is the square of the stretch of the i^th fiber family (=λθ2sin2α0i+λz2cos2α0i); α0i is the fiber angle relative to axial direction in the reference configuration. λ_θ = (r_i + r_o)/ 2ρ_mid, λ_z = l/L are the mean circumferential and axial stretch, and ρ_mid is the unloaded mid-wall radius. The best-fit values of the material parameters and the fiber angle were determined via nonlinear regression of biaxial data from all seven passive protocols. More details on parameter estimation can be found elsewhere (Ramachandra and Humphrey, 2019).

Pulse wave velocity depends on both the geometry and mechanical properties of the arteries and can be well-approximated based on the material stiffness or distensibility of the arterial walls, namely

where ρ is the mass density of the contained fluid (approximately 1,050 kg/m³) and D is the distensibility coefficient (in Pa^–1or kg^–1m⋅s²) determined by the normalized change in arterial diameter from end-diastole to end-systole divided by the change in end-diastolic and end-systolic pressures (Bramwell and Hill, 1922). PWV in the pulmonary artery and its effects on right ventricular hemodynamics can thus be measured non-invasively using the aforementioned variables (Peng et al., 2006; Sanz et al., 2008; Gupta et al., 2018).

For variables with only two categories, differences between hypoxic and normoxic groups in morphological, histological, mechanical, and contractile properties were determined by a two-tailed, unpaired t test. For categorical variables with more than two categories, levels of stretch or pressure were compared using a two-factor analysis of variance (ANOVA). A global test across all levels was performed and then pair-wise comparisons were conducted with post-hoc tests using Bonferroni adjustment. A p < 0.05 level of significance was used, with data reported as mean ± standard error from the mean (SEM).

Viable cells from the main pulmonary arteries of one mouse exposed to hypoxic conditions for six weeks and one age-matched normoxic mouse were analyzed. Following euthanasia, the heart, lungs, and pulmonary arteries were excised with surrounding tissue. The pulmonary arteries were chopped mechanically and placed in 1 mg/ml collagenase (Roche) and 3 U/ml elastase (Worthington). Following cellular dissociation, we barcoded unique mRNA molecules of each cell using our 10× Genomics Chromium platform (3’ v3.1 kit), a droplet-based microfluidic system, and performed reverse transcription, cDNA amplification, fragmentation, adaptor ligation, and sample index PCR according to the manufacturer’s protocol. High sensitivity DNA bioanalyzer traces of cDNA after barcoding and of the final cDNA library were evaluated for quality control. The final cDNA libraries were sequenced on a HiSeq 4000 Illumina platform in our core facility aiming for 150 million reads per library. Raw sequencing reads were demultiplexed based on sample index adaptors, which were added during the last step of cDNA library preparation. Possible adaptor and/or primer contamination were removed using Cutadapt. We processed the RNA-sequence data using the scRNA-seq implementation of STAR (STARsolo), where reads were mapped to the murine reference genome GRCm38 release M22 (GRCm38.p6), collapsed and counted, and summarized to a gene expression matrix. Data were analyzed and visualized using the R packages Seurat (Stuart et al., 2019). Specifically, we clustered the cellular transcriptomes and visualized them in a uniform manifold approximation and projection (UMAP) space to delineate cell types. To identify aberrant gene expression profiles in cellular subpopulations, we established differentially expressed genes between normoxic and hypoxic cells using the non-parametric Wilcoxon rank sum test with Bonferroni adjustment of multiple testing. We used connectomic analyses (R package “Connectome”) (Raredon et al., 2019) to identify conserved patterns of cell-cell cross-talk between the ECs, SMCs and fibroblasts in the pulmonary artery (Raredon et al., 2019; Browaeys et al., 2020). Connectomtic analysis uses an algorithm to visualize intercellular signaling using ligand-receptor mapping. It filters single-cell data to identify the distribution of ligand and receptor z-scores and percent expression to identify cell-type specific communication patterns (edges) with increased statistical confidence. The thickness of the edges corresponds to the z-score (Raredon et al., 2019). Enriched terms analyses of phenotype (using MGI Mammalian Phenotype Level 4 2019 database) and signaling pathways (using WikiPathways 2019 Mouse database) were based on differentially expressed genes when comparing hypoxic with normoxic genetic expression (Chen et al., 2013; Kuleshov et al., 2016). The top results of the gene-set libraries of the collective functions of gene lists generated from the gene list of our single cell sequencing experiments were displayed and plotted as bar charts to show the top enriched terms in the chosen library based on their p-value (plotted as −log₁₀p-value).

Biomechanical phenotyping of normoxic (n = 6) and hypoxic (n = 5) RPAs revealed a general structural stiffening in hypoxia despite modest changes in wall thickness (Figure 2). The mean pressure-diameter findings are shown at the vessel-specific optimal, or in vivo, axial stretch and similarly for the circumferential stress-stretch findings; the axial stress-stretch findings are shown at a common pressure of 15 mmHg. The increased structural stiffness was reflected further by a significant decrease in distensibility and associated increase in the calculated PWV. The modest decrease in wall thickness in hypoxia (inferred during ex vivo testing and confirmed via multiphoton microscopy; Supplementary Figure 1) implicated an increase in material stiffness in hypoxia, which in turn suggested a change in extracellular matrix. Material stiffness and elastically stored energy were computed using the four-fiber family constitutive relation, with associated best-fit material parameters (Supplementary Table 1). Calculated values of biaxial material stiffness (circumferential and axial) were greater for hypoxic than for normoxic vessels in the axial direction alone at diastolic, mean, and systolic pressures (5, 15, 25 mmHg, respectively), as shown in Figure 2. Hypoxia did not change the elastically stored energy at diastolic or mean pressures but did decrease it at systolic pressures (Supplementary Figure 2).

Wall composition and characteristics of intramural collagen are shown in Figure 3 for representative RPAs from normoxic and hypoxic mice. Movat-stained cross-sections of normoxic and hypoxic vessels suggested no major changes in composition (Figures 3A–C). Area fractions for normoxic and hypoxic RPAs were computed for elastin (37 ± 3.5% vs. 41 ± 4.4%, respectively, p = 0.43), collagen (42 ± 8.0% vs. 32 ± 7.8%, p = 0.40), cytoplasm (12% fraction ± 2.9% vs. 17 ± 3.3%, p = 0.27), ground substance (7% fraction ± 0.9% vs. 6 ± 0.8%, p = 0.20), and fibrin (2% fraction ± 0.3 vs. 4 ± 1%, p = 0.14). Pico-sirius red stained sections imaged with polarized light (Figure 3D) revealed, however, that medial collagen was reduced in hypoxic compared to normoxic vessels (p < 0.05). Alignment and orientation of collagen was assessed using multiphoton imaging, with characteristic images shown at 15 mmHg for both groups (Figures 3E–H); additional images are included in Supplementary Figure 3 and Supplementary Material. The mean orientation of collagen fibers (primarily adventitial) was less axially directed for the hypoxic (10.8° ± 5.3°) than the normoxic (2.0° ± 2.1°) vessels (p = 0.20). Collagen fiber alignment was significantly greater for the hypoxic (κ = 14.3 ± 1.9) than the normoxic (κ = 7.0 ± 1.0) vessels (p = 0.02), where the parameter κ and width of distribution are inversely related (i.e., lower κ equals greater width of distribution and vice versa). To the author’s knowledge, these are the first published measurements of collagen fiber alignment in pulmonary arteries of mice. Note, too, that immunohistochemistry delineated changes in these fibrillar collagens by type, I versus III, revealing a decrease in the former and increase in the latter with hypoxia (Supplementary Figure 4).

Pearson correlation analysis revealed that collagen fiber orientations correlated inversely with distensibility (R² = 0.88), but not wall thickness (R² = 0.07), as shown in Figure 3I. Collagen fiber orientation also correlated directly with axial material stiffness (R² = 0.77) but less so with circumferential stiffness (R² = 0.48), as shown in Figure 3J. Multiphoton imaging also revealed an increase in the density of adventitial cells in the hypoxic compared with the normoxic arteries (Figure 3K). Similar results were found in H&E-stained paraffin samples (Supplementary Figure 5), noting that the ratio of the medial : adventitial cross-sectional areas did not change (Supplementary Figure 6). Immunohistochemistry and immunofluorescence were used to confirm that the layers and morphology of cells of the intimal, medial, and adventitial layers were not co-located (Supplementary Figure 5). Recall that paraffin-embedded histological cross-sections are sampled at ∼5 μm increments along the length of the artery, whereas multiphoton imaging is along an entire 500 μm length. Noting that most of the collagen fibers were in the adventitia, enriched terms analysis of phenotype and signaling pathways are shown for fibroblasts in Figure 3L. The bar chart shows the top enriched term search results based on their p-value (plotted as −log₁₀p-value). Differential gene expression in fibroblasts due to hypoxia associated with phenotypic changes in extracellular matrix previously associated with inflammation and endovascular morphology. Pathways associated with the differentially expressed genes observed due to hypoxia included PI3K-Akt-mTOR, TGFβ, and p38-MAPK, as highlighted in Figure 3M. Differentially expressed genes associated with NF-kB, TLR signaling, and IL-3 signaling were also observed. There was a large increase in the number of ribosomal genes expressed in response to hypoxic exposure.

SMCs endow the pulmonary artery with the ability to vasoconstrict and vasodilate, which were assessed as responses to vasoactive substances at a fixed pressure (15 mmHg) and optimized specimen-specific axial stretch (Figure 4), both over time (panels A and B) and as a percent change based on normalized outer diameter at the end of 15 minutes (C and D). The normalized reduction in outer diameter during vasoconstriction (KCl and PE) of control mice (normoxia) was 31 and 35%, respectively, at 15 min, similar to previous findings (Ramachandra and Humphrey, 2019), while hypoxic exposure (10% FiO₂ for 3–6 weeks) resulted in a statistically significant decrease (p < 0.05) in this measure of contractility ex vivo for both KCl (14%) and PE (28%). Multiphoton imaging of the SMCs yet revealed that their mean orientation was more circumferential (θ–90°) in hypoxic (84°–88°) than in normoxic (63°–86°) vessels (panel E). SMC alignment was also greater in hypoxic (κ = 14.12–20.79) than in normoxic (κ = 7.59–9.03) vessels (panel F). Enriched terms analysis of phenotype and signaling pathways are shown for SMC cells (panel G). SMC differential gene expression due to hypoxic exposure associated with phenotypic changes in morphology, matrix, and cell cycle. Pathways associated with the differentially expressed genes observed due to hypoxia include selenoproteins, proteins of myocontractility, and oxidative stress. Differentially expressed genes associated with PI3K-Akt-mTOR and Ras-MAPK intracellular pathways, as highlighted in Figure 4G. There was a large increase in the number of ribosomal genes expressed in response to hypoxic exposure.

Endothelial cell-dependent vaso-dilatation was measured in response to acetylcholine (Ach), with further functional assessment following eNOS inhibition with L-NAME, the latter shown in panel H. The normalized baseline outer diameter in this case (1.0) is the vessel-specific outer diameter following maximal contraction with PE and after 10 minutes exposure to 10 μM Ach. Panel I reveals no significant hypoxia-induced change in EC induced vasodilatation ex vivo. The authors are not aware of previously published data regarding EC function in the ex vivo mouse pulmonary artery. Multiphoton imaging further revealed an increase in the density of intimal cells of hypoxic compared with normoxic vessels (Panel J, images, and panel K, quantitative). Results from H&E-stained paraffin-embedded sections were consistent (Supplementary Figure 5). Enriched terms analysis of phenotype and signaling pathways for ECs (panel L) revealed that differential gene expression due to hypoxic exposure associated with phenotypic changes in focal adhesion, endovascular morphology, thrombosis, neovascularization, and protection against ischemia. Pathways associated with the differentially expressed genes observed due to hypoxia include TGFβ and PI3K-Akt-mTOR signaling and endothelial processes such as thrombosis (prostaglandins and thromboxane), oxidative phosphorylation, and vasoactivity. These changes and associated structural or functional alterations are summarized in Figure 4L. There was a large increase in the number of ribosomal genes expressed in response to hypoxic exposure.

We identified 3,142 resident cells from a segment of the pulmonary artery from one mouse exposed to hypoxic conditions (FiO₂ 10%) for six weeks and one age-matched control mouse exposed to normoxic conditions (FiO₂ 21%). These cells were represented on a uniformed manifold approximation and projection (UMAP) of cells clustered by genetic markers and cellular origin for normoxia versus hypoxia (Figure 5). Identification of resident cell types was based on their expression of marker genes as shown in Figure 5B: endothelial cells (Dkk2, Gja5, Cxcl12, Bmx, Efnb2, Fbln2, Gja4, Dll4, Hey1, Htra1, Cldn5, Cdh5, and Bmp6), smooth muscle cells (Tagln, Acta2, Myh11, Mustn1), and fibroblasts (Tshz2, Vcan, Col1a1, Col1a2, Col14a1, Fbln1, Dcn, Lum, Aspn, Pi16, Serpinf1, Cygb, and Htra3). There was a decreased proportion of SMCs and increased proportion of ECs in the hypoxic sample, as shown in the inset of Figure 5A. Figure 5C shows heat maps for differentially expressed genes in ECs, SMCs, and fibroblasts. There were 331 differentially expressed genes in hypoxic ECs, with higher expression associated with hypoxia-induced proteins for cellular differentiation and proliferation and angiogenesis such as Igfb7, Ifg1r, Igf2r, Myc, Sox5, Lgals1 and lower expression of genes associated with oxidative respiration metabolism such as Atp5mpl, mt-Atp6Atp5e, and Cox7c. Hypoxic ECs also had lower expressions of intercellular junction proteins such as Icam2, Gja5, Cdh11, and Cldn5 with higher expression of cytoskeletal and extracellular proteins such as Fbln5, Actn1, Eln, and Ccn2, and Csgalnact1. There is higher expression of genes associated with endothelial-mesenchymal transitioning, such as Klf2, Fn1, Bmp6, Smad6, Tcf4, and Eng (Medici and Kalluri, 2012; Chen et al., 2016). Finally ECs had a higher expression of chemotactic proteins such Vwf and Hbegf in hypoxia as well as Edn1, which encodes endothelin-1, a potent vaso-constrictor. There were 107 differentially expressed genes in hypoxic SMCs, including a higher expression of genes associated with cytoskeletal rearrangement such as Tgfb3, Vim, Actn4, Tnnt2, Col18a1, Ecm1, and Smtn (van der Loop et al., 1996; Worth et al., 2001; Chen et al., 2016; Pedroza Albert et al., 2020) and cell signaling such as Gne, Pde10a, and Map3k19. Finally, there were 113 differentially expressed genes in hypoxic fibroblasts, with higher expression of genes associated with cellular adhesion and migration such as Cd44, Rock2 and Fn1 and the degradation of extracellular matrix and the ubiquitination-proteosome system such as Smurf2, Tnfaip3, and Ntn4. This is the first single cell RNA sequencing of mouse pulmonary arteries of which the authors are aware, noting that the induced hypoxia resulted in an increase in right ventricular pressure and cavity volume (Supplementary Figure 7).

Differential connectomic analysis identifies changes in cell-type-specific molecular communications that occur among ECs, SMCs, and fibroblasts when comparing hypoxic and normoxic conditions based on previously described protein interactions (Figure 6). These results demonstrate how differentially expressed genes within one cell type may affect gene expression in another cell type, ultimately altering structure and function at the tissue level. Fn1/Vim-Cd44 connections between ECs and fibroblasts in the RPA involves a pathway that has been associated with ECM-cell receptor interactions (Geng et al., 2019) and endothelial-to-mesenchymal transition (Endo-MT) (Goveia et al., 2020). The Tgfb3-Eng connection between SMCs and ECs is associated with EC adhesion and regulation of ECM formation (Cheifetz et al., 1992; Kaartinen et al., 1995). Hbegf-Cd9 signaling seen in the connections between ECs and fibroblasts has previously been associated with preservation of membrane barrier function (Schenk et al., 2013) and enhancing cell viability (Takemura et al., 1999). Enriched terms analysis of phenotype and signaling pathways based on genes identified through differential connectome analysis revealed associations with abnormal lung vasculature, altered collagen deposition, alterations of cell cycle, regulation of actin cytoskeleton, TGF-beta signaling, and alterations in focal adhesion. This captures phenomena and pathways of importance that we observed in this study. These are the first time that these communications between the resident cell types of the mouse pulmonary arteries resulting from hypoxic exposure have been described in this manner.