Label-free metabolic biomarkers for assessing valve interstitial cell calcific progression

Calcific aortic valve disease (CAVD) is the most common form of valve disease where the only available treatment strategy is surgical valve replacement. Technologies for the early detection of CAVD would benefit the development of prevention, mitigation and alternate therapeutic strategies. Two-photon excited fluorescence (TPEF) microscopy is a label-free, non-destructive imaging technique that has been shown to correlate with multiple markers for cellular differentiation and phenotypic changes in cancer and wound healing. Here we show how specific TPEF markers, namely, the optical redox ratio and mitochondrial fractal dimension, correlate with structural, functional and phenotypic changes occurring in the aortic valve interstitial cells (VICs) during osteogenic differentiation. The optical redox ratio, and fractal dimension of mitochondria were assessed and correlated with gene expression and nuclear morphology of VICs. The optical redox ratio decreased for VICs during early osteogenic differentiation and correlated with biological markers for CAVD progression. Fractal dimension correlated with structural and osteogenic markers as well as measures of nuclear morphology. Our study suggests that TPEF imaging markers, specifically the optical redox ratio and mitochondrial fractal dimension, can be potentially used as a tool for assessing early CAVD progression in vitro.

www.nature.com/scientificreports www.nature.com/scientificreports/ to assess the mitochondrial organization by determining mitochondrial fractal dimension (FD) 18,20 . Previously, TPEF was used to monitor the progression of osteogenic differentiation in human mesenchymal stem cells (MSCs) 18,21,22 and ORR and mitochondrial FD was shown to correlate with the osteogenic differentiation of human MSCs 18 . Recently, TPEF imaging was utilized in vitro and on ex vivo tissue explants to detect mineralization, a key hallmark of CAVD 23 . In the context of VICs, we have previously shown that ORR correlated with cellular proliferative potential when VICs were cultured under different media conditions 24 . We have also previously reported that an increase in pathological stretch reduced the ORR in VICs, suggesting that the involved signaling pathways and VIC pathological function are closely linked to the overall cellular metabolic state 14,19,24 . However, TPEF imaging metrics -ORR and mitochondrial clustering have not yet been shown to predict or correlate with the pathological changes in VICs during early CAVD progression.
The objective of this study was thus to investigate the potential of ORR and mitochondrial organization as label-free markers for tracking early CAVD progression. We seeded porcine aortic VICs as monolayers in quiescent versus osteogenic media on two dimensional soft or stiff substrates. We examined these samples using TPEF microscopy to quantify ORR and mitochondrial FD and simultaneously characterized the CAVD progression in our model using traditional end-point biomarkers, such as calcific nodule quantification, gene expression, cell proliferation and apoptosis. We then correlated ORR and FD with the aforementioned VIC structural and biological metrics. Our results showed that TPEF metrics correlated with the early markers of CAVD progression and thus suggest that TPEF microscopy can be utilized as a label-free non-destructive tool for assessing CAVD progression in vitro.

Results
Osteogenic cultures on compliant and stiff substrates developed calcific nodules over time. Progression of nodule size and number has been used as a classic marker for CAVD progression in vitro 2,9 . In our study, VICs cultured under quiescent conditions showed little to no calcium-rich nodules (Fig. 1). VICs under osteogenic conditions developed nodules over time (Fig. 1a), suggesting robust disease progression in our model. Positive Alizarin Red S (ARS) staining on day 28 confirmed that the nodules present in the osteogenic cultures by day 28 were calcific in nature (Fig. 1b). The number of nodules per sample increased with time (Fig. 1c). The number of nodules in compliant substrate osteogenic cultures at day 28 (p = 0.0356) were higher than the number of nodules in quiescent cultures and compliant substrate osteogenic culture at day 1. The normalized area of nodules increased with time (Fig. 1d), with the area of nodules in compliant substrate osteogenic cultures at day 21 (p = 0.0316) being higher than the area of nodules in quiescent cultures and compliant substrate osteogenic culture at day 1.
Optical Redox Ratio (ORR) decreased at day 14 and then increased again at day 28. Label-free TPEF microscopy was performed on live VICs to generate ORR maps (Fig. 2a). Average ORR decreased over time till day 14 and then increased again by day 28, resulting in a lower ORR at day 14 (p < 0.0001) compared to all other time points (Fig. 2b). Overall, ORR for osteogenic cultures was significantly lower than quiescent cultures (p = 0.0072). VICs cultured on compliant substrates showed more pronounced changes in ORR as compared to those cultured on stiff substrates. Within samples cultured on compliant substrates, ORR observed in osteogenic VICs at day 14 was significantly lower than that of VICs at day 1 (p < 0.0001), day 3 (p < 0.001), day 21 (p < 0.005) and day 28 (p < 0.001). Within samples cultured on stiff substrates, ORR osteogenic VICs at day 14 was significantly lower than that observed at day 1 (p = 0.0109), and increased significantly by day 28 (p = 0.0453). Fractal dimension (FD), a metric inversely proportional to mitochondrial clustering, was assessed from the NADH images (Fig. 2c). FD for osteogenic VICs on compliant (p = 0.002) and stiff (p = 0.0267) substrates was significantly higher than that of quiescent VICs on stiff substrate on day 28.
ORR correlated with VIC osteogenic markers and FD correlated with VIC structural markers. To assess the utility of TPEF markers for tracking CAVD progression in vitro, we correlated ORR (Fig. 3) and FD (Fig. 4) with fold change in gene expression of structural, functional and phenotypic markers using Spearman's rank correlation coefficients. These markers included the genes ACTA2, RHOA, TGFβR1, OPN, OCN and RUNX2. Correlations were performed with respect to early (days 1-7) and late (days 14-28) culture time points, media (quiescent vs osteogenic), substrate (compliant vs. stiff) and their combinations. While ORR correlated mostly with phenotypic and functional markers TGFβR1, OPN, OCN and RUNX2 (Table 1a), FD also correlated with structural markers ACTA2 and RHOA (Table 1b). TGFβR1 codes for transforming growth factor beta (TGFβ) receptor 1, which is upregulated in cells on stiffer matrices and has been shown to steer nodule formation on stiffer matrices via apoptosis instead of osteogenesis 9,25,26 . For all time points and conditions, ORR correlated with TGFβR1 (Sp ρ = 0.5411, p < 0.0001) and these correlations were found to be stronger in stiffer matrices especially in the osteogenic conditions and at later time points (Table 1a). Osteopontin (OPN), a known early biomarker for osteogenic differentiation in VICs 27 , correlated with ORR (Sp ρ = 0.4494, p < 0.0001) for all conditions and time points but showed the strongest correlations at in early osteogenic VICs (Sp ρ = 0.8316, p < 0.0001) (Table 1a). Additionally, RUNX2 is a key transcription factor shown to promote CAVD via multiple effectors 28,29 and its expression also induces expression of other osteogenic markers like (OPN) and osteocalcin (OCN) 2 . ORR also correlated with TGFβR1 (Sp ρ = 0.609, p = 0.0044), OPN (Sp ρ = 0.8316, p < 0.0001), OCN (Sp ρ = 0.8203, p = 0.0003) and RUNX2 (Sp ρ = 0.5398, p = 0.014) gene expression in early osteogenic VICs. Though RUNX2 had stronger correlation with ORR on stiffer substrates while not showing any correlations on compliant substrates, OCN correlated with ORR much stronger on compliant substrates especially in the osteogenic conditions (Table 1a). Finally, the structural markers ACTA2 and RHOA play important roles in cellular organization. ACTA2 encodes alpha smooth muscle actin (αSMA) responsible for generation of stress fibers, and RHOA encodes a molecular switch involved in actin re-organization during nodule formation 30,31 . TGFβ1 is known to www.nature.com/scientificreports www.nature.com/scientificreports/ also induce αSMA expression and myofibroblastic differentiation in VICs, which increases with an increase in substrate stiffness 32 (Table 1b). FD showed stronger correlations with RHOA in early timepoints especially in early osteogenic VICs (Sp ρ = 0.5218, p = 0.0183). FD showed strong correlation with ACTA2 in later time points especially in stiffer matrices (Sp ρ = 0.5588, p = 0.0197), with FD correlating strongly and with TGFβR1 (Sp ρ = 0.502, p = 0.0173) in early VICs especially on stiffer matrices (Table 1b).
Nuclear area and aspect ratio correlated with fractal dimension. We have previously demonstrated that VIC shape influences its functional characteristics such as contractility, metabolism and proliferation 19,33 . We have also shown that VIC aspect ratio under normal and osteogenic culture conditions correlates with nuclear aspect ratio 33 . In the current study, ORR maps were utilized to assess nuclear area ( Fig. 5a) and nuclear aspect ratio (AR) (Fig. 5d) and then correlated with ORR and FD. Nuclear area was significantly lower at days 3, 7, 14, 21 and 28 as compared to day 1 (p < 0.001). Nuclear area was significantly higher for osteogenic VICs as compared to quiescent VICs on both compliant substrates (p = 0.0008) and stiff substrates (p = 0.0017) at day 14, which also interestingly coincided with the lowest ORR. Overall, nuclear area was significantly higher for osteogenic VICs (p = 0.0031) compared to the quiescent VICs, especially on compliant substrates (p = 0.0301). AR decreased significantly at day 3, as compared to day 1 (p = 0.0364) and then increased again at day 7 (p < 0.0001) and day 14 (p < 0.0001). AR decreased at day 21 (p = 0.0062) and day 28 (p = 0.0235), as compared to day 14. At the 14-day time point, osteogenic VICs on complaint substrates had significantly lower AR as compared to quiescent VICs on complaint substrate (p = 0.05). When correlating nuclear area with ORR and FD, we observed that ORR did not correlate with nuclear area (Fig. 5b), but FD positively correlated with nuclear area (R = 0.3470, p < 0.0001) www.nature.com/scientificreports www.nature.com/scientificreports/ (Fig. 5c). The nuclear AR was then correlated with ORR and FD. ORR did not correlate with nuclear AR (Fig. 5e), however, FD negatively correlated with nuclear AR (R = −0.3047, p = 0.0002) (Fig. 5f).
Nodules in osteogenic cultures were apoptotic and proliferative in nature. VICs have been known to undergo apoptosis and proliferation during CAVD progression. Additionally, apoptosis is a known marker for dystrophic calcification 34 . Caspase staining revealed the presence of apoptotic cells in the nodules (Fig. 6a). Mean apoptosis index of the nodules was quantified (Fig. 6b) and was found to be significantly higher in the osteogenic media conditions compared to their corresponding quiescent media conditions (p < 0.05) and correlated with RHOA gene expression (R = 0.4049, p = 0.0497). Ki67 staining revealed that the cells were proliferative within, and outside, the nodules (Fig. 7a). Proliferation index for cells outside of the nodules (Fig. 7b) and those within the nodules (Fig. 7c) was quantified and it was seen that the proliferation index was significantly higher in the osteogenic media conditions, compared to their corresponding quiescent media conditions (p < 0.05). The proliferation index for cells outside of the nodules correlated with OCN gene expression

Discussion
Currently, there are few tools available to non-invasively monitor CAVD progression in vitro, thus somewhat limiting the efficient development of therapies for this disease. In this study, we hypothesized that non-invasive TPEF metrics would correlate with the early phenotypic changes of VICs during CAVD progression. This hypothesis was tested using an in vitro CAVD model via the modulation of media and substrate stiffness in two-dimensional VIC cultures. The presence of calcific nodules, apoptosis and proliferation of cells within nodules, and gene expression were used to assess functional changes in VICs. Nuclear morphology was used to describe the structural properties of VICs. TPEF metrics of ORR and FD were then correlated with these structural and functional markers.
ORR and FD were previously shown to correlate with the osteogenic differentiation of mesenchymal stem cells 18 . In assessing the ORR of VICs cultured under quiescent and osteogenic conditions, we show for the first time that ORR changed over time during early CAVD progression. Specifically, the ORR was significantly lower at day 14 and then increased again at day 21. This trend was similar to that observed during osteogenic differentiation of mesenchymal stem cells in vitro 18 . In our model, more pronounced changes were observed in the ORR compared to mitochondrial reorganization as measured by fractal dimension (FD). It has been previously demonstrated that changes in ORR can occur before changes in mitochondrial organization within the cell 18 , and our results support this concept as well. Additionally, we also reported a correlation between osteogenic gene expression and ORR and FD as seen in prior studies 18 , suggesting the possibility of using TPEF metrics to predict the CAVD disease process.
Digging deeper into our results, we observed that TGFβR1 expression significantly correlated with ORR, which was expected given that TGFβ1 signaling has a major role in inducing disease during valve calcification. RUNX2 expression is known to predict the early osteogenic lineage of the cell 25 , and thus correlated well with ORR during early timepoints. Additionally, RUNX2 and OCN correlating with ORR based on stiffness, further supported the mechanosensitive nature of RUNX2 signaling in VICs. RUNX2 and OCN also correlated with VIC proliferation index, suggesting that osteogenic cells tended to be more proliferative. In the context of the above results, it is also important to note that while our quiescent media and osteogenic media formulations contain differing amounts of FBS, previous reports have rigorously characterized these specific media formulations to www.nature.com/scientificreports www.nature.com/scientificreports/ maintain quiescent and osteogenic cells, respectively 15,24,25,[33][34][35][36][37][38][39] . Our results also agree with those accounted previously that demonstrated that quiescent VICs are less proliferative than VICs in 10% FBS containing osteogenic media 25,36,[38][39][40] . Future studies will focus on the effect of individual media formulations and FBS concentrations on TPEF metrics and cell proliferation.
OPN expression has been known to be an early biomarker for CAVD and has been used as a serum marker to study the progression of the disease when it is asymptomatic 27,41 . In our study, OPN and OCN correlated with the ORR, specifically for the early timepoints. While TGFβR1, RUNX2, OPN, and OCN gene expression all were highly correlated with ORR for early timepoints, OPN expression showed the strongest correlation coefficient. OPN did not show any correlations to ORR in late time points, but correlated well with osteogenic cultures under both stiff and compliant cultures at early timepoints. Taken as a whole, our results suggest that the correlation of ORR with OPN can potentially be exploited as a predictor of early metabolic changes during CAVD progression and serve as a biomarker for early disease 27,41 .
Alpha SMA organizes into stress fibers and affects stiffness, force generation, and cell shape 30 . We have previously shown that elongated cells had higher expression of proliferation and αSMA expression 19 . Dystrophic calcification has been known to be associated with apoptosis and prolonged αSMA activation in VICs 30 . In our study, ACTA2 did not correlate with ORR but significantly correlated with mitochondrial FD. Stronger correlations were observed between ACTA2 and FD in osteogenic VICs on stiffer substrates, especially at later time points. This suggests that the role played by ACTA2 in stress fiber formation may also affect the regulation of the arrangement of mitochondria in adjusting to the energy demands of the cell under altered stiffness and osteogenesis. Stiffness dependent correlations of TGFβR1 with FD at early timepoints may suggest its role in triggering the aforementioned ACTA2 mediated mitochondrial arrangement 32 . Additionally, the RHOA gene, which codes for a mechanotransductive protein that regulates the formation of nodules 42 , also correlated with nodule apoptotic index. RHOA has been known to be activated and expressed at higher levels in stiffer matrices and promote the upregulation of RUNX2 43 . Like ACTA2, RHOA did not correlate with ORR but significantly correlated with mitochondrial FD. RHOA has been shown to increase ACTA2 expression, ultimately increasing stress fiber formation 30,31 . This sequence of events is in line with our results that showed that RHOA gene expression correlated with FD at earlier time points, whereas ACTA2 gene expression correlated with FD for later time points.  Continued www.nature.com/scientificreports www.nature.com/scientificreports/ Due to the impact of VIC shape on cellular function, changes in nuclear area and aspect ratio were measured to provide context for the differences we observed in mitochondrial FD. These measures significantly correlated with mitochondrial FD, with nuclear area positively correlated with FD, and aspect ratio negatively correlated with FD. This suggests that mitochondrial clustering decreased with an increase in nuclear area and a decrease in nuclear aspect ratio. We observed an increase in nuclear area and FD in osteogenic VICs, suggesting these cells had decreased mitochondrial clustering. Similar to our previous study 33 , the VICs in osteogenic conditions tended to have a lower nuclear aspect ratio, reinforcing the idea that osteogenic cells have reduced mitochondrial clustering. Disorganization of mitochondria has previously been associated with other cardiac ailments 44 , but the specific functional role of mitochondrial clustering in valve pathophysiology is yet to be elucidated. It is also not fully understood if the increased area and decreased AR of nuclei in osteogenic VICs is a cause or effect of the disorganization of mitochondria.
Overall, the correlations we observed between nuclear morphology, gene expression and the TPEF metrics of ORR and FD, suggest that VIC phenotype and shape have an effect on the metabolic state of the cell and its mitochondrial arrangements. This study showed that TPEF microscopy and its derived metrics can be used to assess pathophysiological changes in the structure, function and phenotype of VICs undergoing osteogenic differentiation, suggesting that TPEF imaging markers can serve as a label-free, non-destructive tool to study the cellular changes occurring during CAVD progression in vitro. Future studies will dig deeper into the exact metabolic changes that occur during CAVD progression and how they relate to observed alteration in measured TPEF metrics, and whether any of these techniques can be translated in vivo.

Materials and Methods
Valve interstitial cell isolation and seeding. Standard cell culture procedures were followed as published earlier 19,33 . Briefly, porcine hearts were acquired from a local abattoir and aortic valves were utilized to aseptically isolate the valve interstitial cells (VICs). Human fibronectin (Corning, Corning, NY) at a concentration of 100 μg/ml was coated on stiff or compliant PDMS-coated coverslips (Supplementary Table S1), before seeding the cells 33 . VICs were conditioned in quiescent media (Dulbecco's Modified Eagle Medium, with 50 U/ ml penicillin, 50 U/ml streptomycin, and 10 mM HEPES, supplemented with 2% fetal bovine serum (FBS), 10 ng/ ml FGF-2 and 50 ng/ml insulin) 36 for 14 days before seeding onto the coverslips. These quiescent VICs from passages three to seven were then seeded at 500,000 cells/coverslip and cultured in either quiescent or pro-osteogenic (Dulbecco's Modified Eagle Medium, with 50 U/ml penicillin, 50 U/ml streptomycin, and 10 mM HEPES, supplemented with 10% FBS, 10 mM β-glycerophosphate, 100 nM dexamethasone, and 50 µg/ml ascorbic acid) conditions 33 for 1, 3, 7, 14, 21, and 28 days before analyses were performed as outlined below.
Assessment and quantification of calcific nodules. Three to four fields of view (FOV) from five to six biological replicates per time point per condition were imaged using brightfield phase contrast microscopy before they were utilized for other assays at the 1, 3, 7, 14, 21, and 28-day time points. An image editing software, Fiji Image J (National Institutes of Health) was used to manually count and measure the number and area of nodules per field of view. At the 28-day time-point, the samples were fixed in 4% paraformaldehyde (Electron Microscopy Sciences, Hatfield, PA) with 1% Triton X-100 (Electron Microscopy Sciences) and incubated with 40 mM Alizarin Red S (ARS; Sigma-Aldrich, St. Louis, MO) solution for 30 minutes with gentle shaking 45 . The samples were then washed in deionized water and brightfield images were acquired using a Nikon Eclipse (Ci-S) microscope in conjunction with NIS Elements software. www.nature.com/scientificreports www.nature.com/scientificreports/ (Middleton, Wisconsin) inverted for cell culture imaging and equipped with an ultrafast Ti:Sapphire tunable laser source (Mai Tai HP; Spectra-Physics, Santa Clara CA) was used to acquire images via a water immersion objective (20×; NA = 1; 2.5x digital zoom) (Olympus, Tokyo, Japan) and non-descanned GaAsP photomultiplier tubes (PMT) (H10770PB-40, Hamamatsu, Shizuoka, Japan). NADH fluorescence was captured with a 460 (±20) nm band pass filter at 755 nm excitation, and FAD fluorescence with a 525 (±25) nm band pass filter at 860 nm excitation. NADH and FAD autofluorescence intensities were normalized by PMT gain and laser power, with PMT gain normalized to fluorescein concentrations (0.1 µM to 20 µM in Tris Buffer of pH 8) as in previous studies 18,46 . PMT gain and laser power were kept constant and laser power was read after each imaging session 18,19 . A custom MATLAB code was used to generate redox images by pixel-wise computation of fluorescence intensities as FAD/ (NADH + FAD) after normalization. Optical redox ratios (ORR) were calculated using averaged and normalized field intensities of FAD and NADH images as FAD/(NADH + FAD). NADH fluorescence intensity maps were utilized to assess the spatial arrangement of mitochondria using a custom MATLAB code. Specifically, a modified blanket method was utilized to obtain the mitochondrial fractal dimension, a metric defining the closeness of mitochondrial arrangement, on a pixel-by-pixel basis 20 . At least two to four fields of view were taken from five to six separate samples per time point, and per treatment group 18,19 . Quantitative real-time polymerase chain reaction and correlation analysis. Quantitative real-time polymerase chain reaction (qRT-PCR) was used to quantify gene expression of ACTA2, RHOA, TGFβR1, OPN, OCN, and RUNX2 using previously published methods 33 . Briefly, a QIAGEN RNeasy Plus assay kit (Qiagen, Hilden, Germany) was utilized to isolate and purify total mRNA from VIC samples. The mRNA was then reverse transcribed to cDNA using iScript RT supermix (Bio-Rad, Hercules, CA). QRT-PCR was performed using the SsoAdvanced Universal SYBR Green supermix (Bio-Rad) in a CFX96 real-time system (BIO-RAD) to quantify the relative expression of genes (Table 2) with respect to the house keeping gene 18S (ΔCt). All primers were obtained from ThermoFisher Scientific (Waltham, MA). The average of delta Ct for all 4 groups for day 1 was then used to normalize individual delta Ct values for all groups at all timepoints to obtain the ΔΔCt. Fold change www.nature.com/scientificreports www.nature.com/scientificreports/ was then calculated as 2 −ΔΔCt and graphed. Three to seven samples per time point and per treatment group were utilized for the analysis. The fold change of each sample was correlated with its corresponding ORR and FD.

Two-photon excited fluorescence (TPEF) imaging and analysis.
Nuclear morphology analysis. ORR maps were additionally used to quantify nuclear aspect ratio and area. Five to eight cells per image field were analyzed using Fiji Image J using previously published methods 33 . Average nuclear area and aspect ratio per treatment group per condition were then calculated and correlated to the average ORR and FD for that treatment group.
Apoptosis and proliferation analysis. Separate samples were fixed using 4% paraformaldehyde (Electron Microscopy Sciences) with 1% Triton X-100 (Electron Microscopy Sciences), blocked with 10% goat serum (Life Technologies, Carlsbad, CA), and were immunolabeled with antibodies against either cleaved caspase 3 (1:400; Cell Signaling Technology, Danvers, MA) to assess apoptosis, or Ki67 (1: 100; Abcam, Cambridge, MA) to assess proliferation. Samples were then fluorescently labelled with appropriate Alexa Fluor-594 conjugated (1:200; life Technologies) secondary antibody and 4′,6-diamidino-2-phenylindole (DAPI, 1:200; Life Technologies) to stain the nuclei. The samples were mounted onto glass coverslips using Prolong Gold antifade agent (Life Technologies) and imaged using a Nikon epifluorescence microscope in conjunction with NIS Elements software. Technical controls were utilized to assess specificity of the secondary antibody (Supplementary figure S4). (Fiji Image J software was utilized to assess the area of positive caspase expression in the nodules using a constant threshold value. An apoptosis index of the nodules was calculated by normalizing the caspase positive area by the total area of nodules per sample. Fiji Image J software was utilized to assess the number of cells with positive Ki67 expression per field of view. The Ki67 positive cells within the nodules were analyzed separately from the outside the nodules. Proliferation index of the nodules was calculated by normalizing the Ki67 positive cells by the total area of nodules per sample. The proliferation index of the cells was calculated by dividing the total number of Ki67 positive cells by the total cell coverage area. www.nature.com/scientificreports www.nature.com/scientificreports/ Statistical analysis. All data was graphically presented as mean ± standard error of means (SEM). Normally distributed data was analyzed using three-way analysis of variance (ANOVA) and Tukey's post-hoc tests. For the analysis of ORR, FD, nuclear area and aspect ratio, each FOV was nested within the sample and each sample was nested within the stiffness and media condition. Both sample and FOV were considered a random variable. ANOVA on ranks with Dunn's multiple comparisons test was used for non-parametric data. A comparison was considered significant if its p-value was less than 0.05. Pearson's correlation coefficient was used for the correlation analysis. Significant correlations were defined based on a null hypothesis of uncorrelated data (R = 0). If the data was non-normal, Spearman's rank correlation coefficient was used. All the statistical and correlation analyses was performed using JMP and graphing was performed using SigmaPlot (Systat Software Inc) software.  www.nature.com/scientificreports www.nature.com/scientificreports/ www.nature.com/scientificreports www.nature.com/scientificreports/