In-vivo correlations between skin metabolic oscillations and vasomotion in wild-type mice and in a model of oxidative stress

Arterioles in the cutaneous microcirculation frequently display an oscillatory phenomenon defined vasomotion, consistent with periodic diameter variations in the micro-vessels associated with particular physiological or abnormal conditions. The cellular mechanisms underlying vasomotion and its physiological role have not been completely elucidated. Various mechanisms were demonstrated, based on cell Ca2+ oscillations determined by the activity of channels in the plasma membrane or sarcoplasmic reticulum of vascular cells. However, the possible engagement in vasomotion of cell metabolic oscillations of mitochondrial or glycolytic origin has been poorly explored. Metabolic oscillations associated with the production of ATP energy were previously described in cells, while limited studies have investigated these fluctuations in-vivo. Here, we characterised a low-frequency metabolic oscillator (MO-1) in skin from live wild-type and Nrf2−/− mice, by combination of fluorescence spectroscopy and wavelet transform processing technique. Furthermore, the relationships between metabolic and microvascular oscillators were examined during phenylephrine-induced vasoconstriction. We found a significant interaction between MO-1 and the endothelial EDHF vasomotor mechanism that was reduced in the presence of oxidative stress (Nrf2−/− mice). Our findings suggest indirectly that metabolic oscillations may be involved in the mechanisms underlying endothelium-mediated skin vasomotion, which might be altered in the presence of metabolic disturbance.

The term vasomotion indicates rhythmic oscillations of blood vessels diameter responsible for changes in the vascular tone and blood perfusion to tissue 1,2 . The physiological role of these oscillations is unclear, with contrasting studies supporting the development of vasomotion in both diseased and healthy conditions 2 . The origin of the oscillations is mainly determined by the local contraction and relaxation of vascular smooth muscle cells (VSMCs), which act as synchronised pacemakers in the vessels wall [1][2][3] . Although in several arteries (e.g. rat aorta) the oscillations may depend exclusively on the activity of VSMCs 1 , in other vascular beds (i.e. rat mesenteric artery) the most internal layer of blood vessels (endothelium) plays an important modulatory role during arterial motion 4,5 . Vasomotion is more prominent in small resistance micro-vessels 2 such as skin arterioles that represent an optimal sample for the in-vivo examination of this phenomenon reflecting the general health conditions of the cardiovascular system 6 .
For the explanation of the cellular mechanisms causing the oscillations of vessels diameter observed during vasomotion, three different kinds of cellular oscillators have been proposed as possible drivers of this phenomenon: cytosolic oscillator, membrane oscillator and metabolic oscillator 1,2  between cytosolic and membrane oscillators to enhance the VSMCs synchronisation responsible for vasomotion has been largely explored and recognised 1,2,[9][10][11][12] . They act inducing Ca 2+ fluctuations in VSMCs, respectively through the oscillatory release of Ca 2+ from the sarcoplasmic reticulum 1,10,11 , and via the opening activity of voltage-dependent Ca 2+ channels and large-conductance K + channels of the plasma membrane 1,2,12 . Furthermore, Ca 2+ and membrane potential oscillations have been described in endothelial cells (ECs), suggesting that in some vascular beds (i.e. rat mesenteric artery) in response to specific stimuli (e.g. α-adrenergic stimulation) a primary oscillation from ECs and the interaction with VSMCs might be essential for vasomotion [13][14][15] . A third hypothesised driving mechanism is the presence of a metabolic oscillator represented by oscillations in the activity of the glycolytic enzyme phosphofructokinase (PFK) that cause fluctuations in glycolysis, ATP levels and in the activity of plasma membrane ion channels and membrane potential 1 . This mechanism has never been considered important because of limited experimental evidence for the involvement of PFK in glycolytic and Ca 2+ fluctuations 1,2 . However, recent studies have strongly suggested the participation of PFK in reactions responsible for glycolysis and Ca 2+ oscillations [16][17][18] , with a possible driving effect exerted in-vivo by external glucose availability and uptake [18][19][20] . Furthermore, metabolic mitochondrial oscillations have also been reported, mainly associated with fluctuations of the mitochondrial membrane potential ΔΨm that may be mediated by calcium [21][22][23] . The external oxygen required for sustaining the electron transport chain has been proposed as the driving force for the mitochondrial oscillator 18 . All this evidence suggests that the energetic and oxygen tissue requirements might stimulate vasomotion phenomena associated with Ca 2+ oscillations of metabolic origin, which may involve the intercellular communication between heterogeneous groups of cells (e.g. ECs, VSMCs and cells of the tissue surrounding blood vessels).
Metabolic oscillations have been studied in different cell types, i.e. yeast Saccharomyces cerevisiae 19,21 , cardiac myocytes 24 , β-cells 23 , and can be detected indirectly by measuring the fluctuations in the autofluorescence of the intermediate product of energy metabolism NAD(P)H (nicotamide adenine dinucleotide). However, there are limited studies describing metabolic oscillations in live tissue. Mitochondrial fast oscillations have been reported in-vivo in rat, and detected by intravital two-photon (2P) microscopy measurements of NAD(P)H 25 . Nevertheless, to our knowledge, the slow metabolic NAD(P)H oscillations described in many in-vitro or ex-vivo cell studies 16,23 have not been investigated in-vivo.
In the present study, we have utilised skin tissue from live mice to simultaneously measure NAD(P)H and microvascular blood flow respectively by laser fluorescence spectroscopy (LFS) and LDF, under physiological conditions and during α-adrenergic stimulation with the vasoconstrictor drug phenylephrine (PE). Data were processed by continuous wavelet transform (CWT) spectral analysis to characterise heterogeneous oscillations of NAD(P)H and LDF signals. We found low-frequency oscillations of NAD(P)H autofluorescence with periods of 2.4 min, 5 min and 10 min. To our knowledge, this is the first study in live skin tissue of slow NAD(P)H oscillations previously described in cells. Furthermore, we found relevant correlations between metabolic and vasomotion oscillations suggesting a link between these processes. Therefore, our results support the idea of vasomotion as a dynamic process in the context of a metabolically active microvascular network involving the cooperation of different cell types and vascular segments 2,15 , and a possible driving mechanism of metabolic origin.

Results
Nrf2 −/− (Nuclear factor-erythroid 2 p45-related factor 2) knockout model. Experiments were performed on 5 wild-type (WT) control and 6 Nrf2 −/− mice. The Nrf2 transcription factor is a key regulator of the cell redox state, involved in the resistance to oxidative stress by activating the cellular antioxidant defence 26 . Nrf2 knockout leads to abnormal mitochondrial activity, altered redox status, increased oxidative stress, and has an impact on cardiovascular dynamics 26,27 . We used the Nrf2 −/− strain to investigate the differences in metabolic and vascular dynamics between normal mice and animals with impaired antioxidant defence that may favour the development of metabolic and cardiovascular disorders. Figure 1 shows the experimental setup (Fig. 1a), and examples of simultaneous LFS and LDF signals (Fig. 1c,d) from the flank of a mouse at baseline (10 min) and during vasoconstriction induced by local administration of PE (10 min). The reason for using PE is that, according to Okazaki et al. 5 , this vasoactive agent is able to stimulate a vasomotion phenomenon mediated by endothelial cells. Thus, the use of PE provided the opportunity to study at the same time the role of the endothelium in vasomotion and its relationship with metabolic oscillations. Because of the small number of animals available in the study, the experiments were repeated in duplicate for each mouse and the average value of each variable between the two experiments was taken as final value for data analysis. Data of all the variables analysed in this study are summarised in Table 1.

General microvascular and metabolic biomarkers.
Results of blood perfusion and oxygen saturation (SO 2 ) are shown in Fig. 2a,b. Blood perfusion decreased to a similar degree in WT (p = 0.046) and Nrf2 −/− (p = 0.007) mice during PE-induced vasoconstriction, however, a stronger p-value was found in knockout mice. SO 2 , measured by reflectance spectroscopy (RS), was similar between groups and decreased slightly following vasoconstriction without statistical significance.
LFS skin irradiation by UV light stimulates the autofluorescence emission of multiple fluorophores, including NAD(P)H (490 nm), flavin adenine dinucleotide (FAD + , 550 nm) and elastin (450 nm) (Fig. 1b). Although the NAD(P)H central emission peak represents the highest contribution to the spectrum, we took advantage of the overlapping contributions from FAD + and elastin, respectively to estimate the redox ratio (RR) index and to normalise NAD(P)H fluorescence: normalised amplitude amplitude RR index is a measure of the mitochondrial redox state, reflecting the balance between the reduced NAD(P) H and the oxidised FAD + coenzymes. Comparing this index between WT and Nrf2 −/− mice should provide discriminatory data, confirming that the Nrf2 −/− phenotype is characterised by altered mitochondrial function due to impaired antioxidant defence and higher oxidative stress. NAD(P)H autofluorescence was normalised by elastin's fluorescence peak to reduce blood volume artefacts (see more details in the Supplementary Information). During each experiment, 20 UV spectra were collected in 20 min at a rate of a spectrum every minute, 10 at baseline and 10 during PE administration. RR index and NAD(P)H normalised values were estimated from each spectrum and used to reconstruct a time series (Fig. 1c) for applying the CWT spectral analysis. Figure 2c,d describe the average trends of NAD(P)H normalised and RR index signals. NAD(P)H normalised did not show significant statistical differences between groups, while RR index was significantly lower in Nrf2 −/− mice (p = 0.016) confirming that knockout animals were affected by altered redox state and higher oxidative stress.

Reproduction of the LDF oscillators reported in the literature. Vasomotion oscillations were
assessed by combining LDF and CWT techniques 8,28,29 . We found oscillatory frequency intervals related to tissue/cell activity similar to those reported in human studies: myogenic (50-150 × 10 −3 Hz), neurogenic (20-50 × 10 −3 Hz), endothelial NO-dependent (9-20 × 10 −3 Hz), endothelial NO-independent (5-9 × 10 −3 Hz) (Fig. 3a). In contrast, cardiac (1350-5000 × 10 −3 Hz) and respiratory (150-1350 × 10 −3 Hz) components showed higher frequency ranges (Fig. 3a), confirming the faster heart and breathing rates expected for mice 30,31 . Although the magnitude of the oscillations could be reduced due to the effects of anaesthesia 32 , we used a light isoflurane anaesthesia to avoid major systemic effects on the circulation. Additionally, all mice were scanned under the same conditions thus making it possible to compare differences between groups.
CWT allows the identification of two oscillators reflecting ECs activity. While the biological link between the 9-20 × 10 −3 Hz oscillator and the endothelial NO mechanism has been demonstrated 33 , the origin of the 5-9 × 10 −3 Hz oscillator (endothelial NO-independent) is ambiguous. However, our results and literature reports 28,29 suggest indirectly an association with the endothelial-derived hyperpolarizing factor (EDHF) vasodilation mechanism, thus here we will refer to the endothelial NO-independent component as EDHF. The reasons supporting this hypothesis are explained in the Supplementary Information.

Parameters
WT (Mean ± SD) Nrf2 −/− (Mean ± SD)  Table 1.  Figure 3b shows an example of CWT scalogram and the corresponding time-averaged spectrum from NAD(P)H signal. High spectral energy was found in the low-frequency ranges, thus slow oscillators dominate the reconstructed tracing. Although time resolution is scarce in the low frequencies, we clearly identified three heterogeneous oscillatory intervals defined metabolic oscillator-1 (MO-1) (5-9 × 10 −3 Hz), metabolic oscillator-2 (MO-2) (2.5-5 × 10 −3 Hz) and metabolic oscillator-3 (MO-3) (1.5-2.5 × 10 −3 Hz). However, we will discuss only results of MO-1 because MO-2 and MO-3 frequency ranges fall largely at the edges of the cone of influence region in the CWT scalogram, where data might not be reliable 34 (see caption of Fig. 3 for more details). Figure 4 describes the quantitative results of the relative spectral energy e i and frequency f extracted from the time-averaged wavelet peak of each LDF oscillator. The relative e i , to allow comparisons between groups, was obtained by normalising the absolute energy extracted from the peak of each oscillator (E i , area under the curve) by the total absolute energy E tot of the wavelet spectrum and the number of frequencies in the interval of the oscillator.

LDF oscillators' changes during vasomotion and comparison between groups.
The spectral e i and frequency of the EDHF endothelial oscillator were significantly greater in WT mice at baseline (energy p = 0.03, frequency p = 0.005). The EDHF energy and frequency changed significantly in both groups during PE stimulation following opposite trends. The values increased in Nrf2 −/− mice (energy p = 0.01, frequency p = 0.02) and decreased in controls (energy p = 0.008, frequency p = 0.04). The endothelial NO oscillator displayed the same trends observed for the EDHF oscillation. However, we found statistical significance only for the decrease of energy during phenylephrine delivery in WT mice (p = 0.008). The neurogenic and myogenic components showed patterns similar to those of the endothelial oscillations. The energy was reduced in normal mice after PE administration, with a good p-value for the neurogenic oscillator (p = 0.04) and a p-value close to statistical significance for the myogenic component (p = 0.08). We evaluated also the relative maximal amplitude a i of the wavelet peaks, estimated by normalising the absolute amplitude A i of each oscillator by the total amplitude A tot of the wavelet spectrum and the number of frequencies in the interval of the oscillator. The data of a i are summarised in Table 1 and displayed the same patterns found for e i but with lower statistical significance.
These results outline opposite microvascular reactivity to PE in WT and Nrf2 −/− mice, especially for the endothelial components that are key factors determining microvascular function. This suggests that oxidative stress may affect endothelial function and the overall cardiovascular dynamics in the long term.
In general, the increase of LDF oscillators' spectral energy reflects vasodilation dynamics 28,29,33 , thus the expected oscillatory pattern during mild vasoconstriction (1% PE) should be a decrease of e i , as observed in WT animals. In contrast, Nrf2 −/− mice displayed an increase of LDF oscillators' e i during phenylephrine administration. The reason for this behaviour might be the activation of an endothelium-mediated vasomotion resistance mechanism to attenuate vasoconstriction defined myoendothelial feedback 15 , consistent with the alternation of vasoconstriction and vasodilation induced respectively by simultaneous targeting of VSMCs with PE and EDHF/ NO endothelial signalling. The cause for the activation of this mechanism may be a general major vasoconstriction in micro-vessels affected by oxidative stress. We cannot exclude that also WT mice may activate this mechanism in response to higher vasoconstriction stimuli (PE > 1%).  The results of the spectral analysis outline different dynamics of metabolic oscillators in WT and Nrf2 −/− mice, especially for the RR MO-1. Considering that RR index reflects mitochondrial function, the major discriminatory power of RR MO-1 may indicate abnormal mitochondrial energetic dynamics in Nrf2 −/− mice due to increased amounts of reactive oxygen species (ROS) affecting the electron transport chain.

Wavelet phase coherence (WPCO). The WPCO analysis allows investigating the phase relationship
Cφ(ωk) between oscillations in a specific frequency range of two signals measured simultaneously. We used WPCO to explore the phase interaction between the endothelial EDHF oscillator and the MO-1 of NAD(P) H normalised (Cφ NAD(P)H MO-1/EDHF ) or RR index (Cφ RR MO-1/EDHF ) signals, which fall in the same low-frequency interval (5-9 × 10 −3 Hz). Because our analysis does not include the surrogate data testing described recently by Gruszecki et al. 35 , we cannot make reliable assumptions on the degree of coherence at rest due to bias affecting the low frequencies. However, we can cautiously compare differences between mice phenotypes and changes in response to PE stimulation. Non-significant differences were found for the Cφ NAD(P)H MO-1/EDHF (Fig. 5d). In contrast, baseline Cφ RR MO-1/EDHF was significantly lower in Nrf2 −/− mice (p = 0.005), and changed significantly in both groups during PE administration showing a decrease in WT animals (p = 0.01) and an increase in Nrf2 −/− mice (p = 0.003) (Fig. 5d). These results outline again a major discriminatory power of RR MO-1 compared to NAD(P)H MO-1, suggesting that oxidative stress may impact the interaction/coupling between mitochondrial reactions and the EDHF vasodilation mechanism affecting vascular reactivity to PE.

Relevant correlations.
Correlations in WT mice. Figure 6 displays the correlations between microvascular and metabolic variables observed in WT mice. Due to the small number of animals (n = 5) available in the study, separating baseline and PE data in two distinct subgroups would reduce the statistical power, thus we have pooled all the data in a single correlation analysis. However, for completeness of information the results are presented for the pooled group (baseline + PE, n = 10) and the separate baseline and PE subgroups (n = 5, statistical power < 0. Correlations in Nrf2 −/− mice. Figure 7 displays the most important correlations observed in Nrf2 −/− mice. Also in this case due to the small number of animals (n = 6), the results are presented either as pooled group (baseline + PE, n = 12) and separate baseline and PE subgroups (n = 6). The main correlations found in the pooled group were related to the wavelet frequency f of the endothelial NO oscillator, which was positively correlated with NAD(P)H MO-1 e i (r = 0.60, p = 0.04) and negatively correlated with SO 2 (r = −0.60, p = 0.03). These correlations might reflect the overexpression of the endothelial nitric oxide synthase (eNOS) enzyme described previously in Nrf2 −/− mice 27 , leading to major production of NO probably to balance for dysfunction in other vasodilation pathways preserving vascular function in vessels affected by oxidative stress 27 . The correlations may be consistent with this evidence because both NADPH and molecular oxygen participate in the reaction for NO biosynthesis 36 . Therefore, the increase of endothelial NO frequency correlated with the decrease of SO 2 and the increase of NAD(P)H MO-1 e i might reflect respectively growth of NO production, consumption of molecular oxygen and oxidation of NADPH in NADP + . We could speculate that a sustained NO production may compensate for dysfunction in the EDHF vasodilator activity that was found significantly lower at baseline in mice affected by oxidative stress compared to WT (Fig. 4a). However, further research with a greater number of animals and a more appropriate model for the study of the NO pathway (e.g. eNOS knockout mice) is required to elucidate better these aspects. Indeed, the results of the baseline and PE subgroups showed a significant correlation between the endothelial NO frequency and NAD(P)H MO-1 e i only at baseline (r = 0.70, p = 0.05), while the endothelial NO frequency and the oxygen saturation were significantly correlated only during PE delivery (r = −0.78, p = 0.04), which may indicate that the two correlations are not linked to the same physiological process.

Discussion
In this work, cutaneous vasomotion and cell metabolic oscillations were examined simultaneously from live mice. A vasomotion phenomenon specifically mediated by NO and EDHF endothelial mechanisms 5 was stimulated through administration of low-dose phenylephrine. Comparisons between Nrf2 −/− mice affected by oxidative stress and WT controls revealed an opposite behaviour of vascular oscillations in the two models, indicating an effect of oxidative stress on microvascular reactivity. The most relevant differences were found especially for the endothelial EDHF oscillator, and vasomotion was more prominent in Nrf2 −/− mice probably due to the activation of the myoendothelial feedback 15 vascular modulation in response to oxidative stress.
We characterised a low-frequency metabolic oscillator (MO-1) of NAD(P)H normalised and RR index signals that displayed different dynamics in WT and Nrf2 −/− models, probably reflecting an alteration of the mitochondrial energetic processes due to oxidative stress. Relevant correlations were found in WT mice between metabolic and microvascular oscillators, highly suggesting an involvement of the cellular processes associated with fluctuations of NAD(P)H concentrations (i.e. OXPHOS and glycolysis) in microvascular reactivity. These correlations were absent in Nrf2 −/− mice suggesting an influence of oxidative stress on the interaction/coupling between cell metabolic reactions and microvascular regulatory mechanisms.
Our findings indicate indirectly that cell metabolic oscillators may have an important role in modulating vasomotion in response to specific stimuli at least for oscillatory processes endothelium-mediated. Therefore, this raises interest in the study of cell mitochondrial and glycolytic oscillators to elucidate further the mechanisms driving the Ca 2+ oscillations at the basis of vasomotion. Despite ECs seemed the microcirculation component mostly associated with slow metabolic oscillators, metabolic fluctuations were also correlated with neurogenic and myogenic components indicating a global cooperation of multiple factors during vasomotion with higher/ lower contributions depending on the type of stimulus and vascular bed. In the context of skin microcirculation, primary metabolic Ca 2+ oscillations coming from ECs and the EDHF may represent the key drivers and regulatory factors for vasomotion. Based on our results, vasomotion can be considered an adaptive mechanism directed by the energetic requirements/stimuli in the environment surrounding blood vessels aimed at ensuring an optimal intake of nutrients (i.e. glucose) and oxygen for the production of ATP energy, rather than a process leading exclusively to the spontaneous pacemaker activity of VSMCs. This definition could help addressing better several controversial aspects, i.e. the observation of vasomotion in both healthy and diseased conditions 1,2 , or the primary involvement of ECs or VSMCs in different vascular beds 1,2 . These contrasting observations may be the result of variable energetic requirements and molecular signals in the parenchymal tissue surrounding the micro-vessels, depending on the degree of nutrition/oxygenation, the presence of a particular pathology, and the cell types involved. According to the literature, the factors released by microvascular ECs (NO, EDHF, etc.) can diffuse to the underlying parenchymal tissue and modulate the mitochondrial metabolism, the production of ROS and inflammation 37 . Thus, considering that the autofluorescence signals measured in this work have a prevalent epidermal origin, the correlations found between metabolic oscillations and vasomotion might also reflect the communication between skin cells and the microvascular network to fulfil the cutaneous metabolic needs in response to the α-adrenergic stimulus. A shortcoming of the present research is the impossibility to extend the results on a general basis due to the small-scale design of the study and the restriction of the experiments exclusively to the WT and Nrf2 −/− models. As future perspective to increase the applicability on a broader basis, further studies are required on a larger number of animals and on additional mouse models. In particular, testing a larger number of Nrf2 −/− mice after challenging with high fat diet or after inducing diabetes would be useful to characterise the behaviour of the nonlinear biomarkers in the presence of an extensive oxidative stress load or a metabolic disease, and testing the biomarkers in a knockout model for the eNOS enzyme would be useful to elucidate better the role of NO.
The experiments in the present work were restricted to female mice due to limited availability of animals, thus another point to address in the future is also the extension of the study to male animals. The groups of tested mice were perfectly matched for an heterogeneous age range (36-60 weeks), which allowed making only genotype-related comparisons. However, this does not represent an issue because the main purpose of the current study was exploring the relationship between metabolic and microvascular oscillations in normal physiology and then use an example where metabolic function was disturbed (Nrf2 −/− model). As future perspective, a longitudinal study using more homogeneous aged animals could be also of interest to characterise the patterns of the nonlinear biomarkers for different age ranges.
The methodology used in this work can be easily translated for the concurrent study of metabolic and microvascular dynamics in humans and the exploration of metabolic oscillators as cardiovascular risk factors for diagnostic purposes. The use of LFS to monitor NAD(P)H provided some advantages, i.e. label-free minimally invasive nature of the measurements, convenient combination of LFS and LDF probes for simultaneous recordings, portability of the technology. However, there are limitations leading to the discrete nature of LFS measurements, necessity to reconstruct NAD(P)H normalised and RR index signals to characterise the oscillators, and the inability to trace the specific cellular origin of the oscillations. Furthermore, the technique was restricted to the investigation of low-frequency fluctuations due to discrete slow sampling of the spectra. More advanced technologies for the in-vivo imaging of NAD(P)H, i.e. fluorescence lifetime-based methods [38][39][40] and multiphoton microscopy 25,[41][42][43] , allow a better monitoring and discrimination between glycolytic and mitochondrial NAD(P) H with high spatial-temporal resolution at cell level. Combining these methods with CWT can help studying both slow and fast NAD(P)H oscillators. Nevertheless, these technologies still require further implementation for the application on humans and for providing continuous NAD(P)H recordings, they are expensive and lacking of portability and device design for combining with the LDF probe to study simultaneously metabolic and microvascular oscillators.  Skin preparation and anaesthesia. LFS and LDF techniques are sensitive to light absorption by hair and skin pigmentation. Therefore, a hair-free non-pigmented intact skin is required to obtain reliable measurements. Forty-eight hours prior to performing the experiments, hair from the mice flanks was shaved using an electric shaver, and the residual hair was removed using depilatory cream (Veet, Reckitt-Benckiser). Before collecting measurements, mice were anaesthetised through a standard Boyle's Apparatus to prevent movement artefacts and were laid on a heat mat at 37.0 °C, as described previously by us 44 . A light anaesthesia was maintained by delivering 1.5-2% isoflurane (Abbott Laboratories) in oxygen (1.5 L/minute) through an inhalation nose cone (Fig. 1a).
Iontophoresis. Iontophoresis allows the local delivery of vasoactive drugs in the skin microcirculation without inducing systemic effects. The drug is transferred transdermally by the unidirectional movement of ions in a solution, through a continuous current generated by a reference electrode and a platinum electrode incorporated in a ring-shaped chamber 44,45 . We used an iontophoresis chamber (ION6 probe, Moor Instruments, UK) of 20 mm internal diameter attached to the mouse flank using double-adhesive tape, and a reference electrode placed underside of the animal. The chamber was filled with a 2 ml 1% phenylephrine solution, and LFS and LDF probes were placed at two adjacent sites on the skin area inside the chamber (Fig. 1a). LFS/LDF signals were measured simultaneously for 20 min, 10 min without application of electric current (baseline), and 10 min during delivery of PE by applying a continuous 100 μA anodal current through a controller (MIC2, Moor Instruments, UK) connected to the electrodes.
Laser Fluorescence Spectroscopy. LFS employs a low-power laser (1-5 mW) for the in-vivo excitation of cell endogenous fluorophores, i.e. NAD(P)H, allowing the detection of autofluorescence emission peaks proportional to the tissue concentration of the biomarkers [46][47][48] . We assessed skin NAD(P)H autofluorescence with a single-point LFS probe (LAKK-M, Spe Lazma, Russia) provided with irradiating and detection sensors at a distance of ~1 mm for targeting a tissue volume of ~1 mm 3 . Discrete 10 seconds autofluorescence spectra were sampled from the mouse flank by a 365 nm UV excitation light (1 mW) at a rate of one spectrum every minute during iontophoresis. The spectra were processed by Matlab R2015a (The MathWorks Inc.) to extract the amplitudes of NAD(P)H (490 nm), FAD + (550 nm), and elastin (450 nm) peaks (Fig. 1b). The redox ratio and normalised NAD(P)H autofluorescence were estimated respectively according to equations 1 and 2. For the study of NAD(P)H normalised and RR index oscillations by CWT, 20 min signals were reconstructed from the discrete data points (Fig. 1c) using the piecewise cubic spline interpolation method that we have previously described for reconstructing blood flow signals 49,50 . More details on signals reconstruction and NAD(P)H autofluorescence normalisation are provided in the Supplementary Information.
Laser Doppler Flowmetry and oxygen saturation. LDF allows the indirect measurement of blood perfusion by sensing the Doppler-shift in wavelength generated when a monochromatic and coherent light is backscattered by moving blood cells 7,8,28,33,45 . We measured skin blood flow using a laser Doppler flowmeter (LAKK-M, Spe Lazma, Russia) provided with a probe identical to the LFS probe described above, which was calibrated with a fluoroplastic oscillating disk simulating Brownian motion. A 20 min LDF signal was measured during iontophoresis (Fig. 1d) using a 1064 nm laser with sampling frequency of 20 Hz. The multi-functional probe was also provided with 630 nm and 532 nm lasers, which allowed measuring the percentage skin oxygen saturation (SO 2 ) according to the reflectance spectroscopy principles 51 . The technique is based on the analysis of the signal backscattered from the skin to detect the differential absorption of red and green visible lights by the oxygenated (oxyHb) and deoxygenated (deoxyHb) haemoglobin (Hb). The absorption is proportional to the relative molar concentrations of Hb fractions in the tissue, which are extracted from the backscattered spectrum by application of the Beer-Lambert law 52 : a i i i where ε i (λ) is the molar extinction coefficient of each absorber, C i is the molar concentration of each absorber, and μ a (λ) is the absorption coefficient. The percentage SO 2 is finally estimated as follow 52 : 2 More details and references on how the terms of equations 3 and 4 can be extracted from the backscattered signal are provided in the paper by Dunaev et al. 52 .

Continuous Wavelet Transform.
The CWT method is largely employed to analyse blood perfusion vasomotion dynamics in the time-frequency domain 8,28,29 , providing optimal frequency resolution for low frequencies and good time resolution for high frequencies. We employed this method to investigate heterogeneous oscillations of LDF, NAD(P)H normalised and RR index signals, using a CWT with scaling factor s, time t, Morlet wavelet function ψ and central frequency 1 8,28 : −∞ ∞ g s t s u t s g u du ( , ) 1 ( ) (5) CWT and the extraction of the spectral energy E i 28 , amplitude A i 28 and frequency f 28 from the time-averaged wavelet peaks were performed using Matlab R2015a (The MathWorks Inc.). Figure 3 shows examples of CWT scalograms and corresponding time-averaged spectra from LDF and NAD(P)H normalised signals.
Wavelet phase coherence. The WPCO analysis provides information on the phase relationship Cφ(ωk) between oscillators in the same frequency interval of two signals recorded simultaneously 53 . The analysis returns values between 0 and 1, where Cφ(ωk) ≈ 0 indicates absence of coherence, Cφ(ωk) ≈ 1 complete coherence, and 0 < Cφ(ωk) < 1 partial coherence 53 . We investigated the coherence between MO-1 of metabolic signals and the LDF endothelial EDHF oscillator that are located in the same frequency range. The analysis was performed through the WPCO code provided in the Lancaster University website http://py-biomedical.lancaster.ac.uk by Professor Aneta Stefanovska research group, which estimates Cφ(ωk) according to the principles described by Bandrivskyy et al. 53 and Clemson et al. 54 .
Statistics. Statistical analyses were performed by R-Studio software. The Shapiro-Wilk test revealed normal distributions for all the variables tested in the study. Therefore, the correlations were evaluated by calculating the Pearson's correlation coefficient r, differences between baseline and PE time points were tested by paired t-test, and differences between groups were assessed by unpaired t-test. The correlations were considered relevant for p ≤ 0.05 and values of r > 0.5 or < −0.5. The t-test was considered significant at p < 0.05 with a statistical power of 0.8.