Transcriptomic Analysis of a Diabetic Skin-Humanized Mouse Model Dissects Molecular Pathways Underlying the Delayed Wound Healing Response
Abstract
This article belongs to the Section Molecular Genetics and Genomics.
Full text
genes G C A T T A C G G C A T Article Transcriptomic Analysis of a Diabetic Skin-Humanized Mouse Model Dissects Molecular Pathways Underlying the Delayed Wound Healing Response Carlos León1,2,3,4, Francisco García-García5,6, Sara Llames 2,3,7 , Eva García-Pérez 7, Marta Carretero 2,3,4, María del Carmen Arriba 1,3,4, Joaquín Dopazo 8,9,10,11 , Marcela del Río1,2,3,4, María JoséEscámez 1,2,3,4,* and Lucía Martínez-Santamaría1,2,3,4 Citation: León, C.; García-García, F.; Llames, S.; García-Pérez, E.; Carretero, M.; Arriba, M.d.C.; Dopazo, J.; del Río, M.; Escámez, M.J.; Martínez-Santamaría, L.; et al. Transcriptomic Analysis of a Diabetic Skin-Humanized Mouse Model Dissects Molecular Pathways Underlying the Delayed Wound Healing Response. Genes 2021,12, 47. https://doi.org/10.3390/genes 12010047 Received: 25 November 2020 Accepted: 23 December 2020 Published: 31 December 2020 Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations. Copyright: © 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/). 1 Department of Bioengineering, Carlos III University, Av. de la Universidad, 30, Leganés, 28911 Madrid, Spain; [email protected] (C.L.); [email protected] (M.d.C.A.); [email protected] (M.d.R.); [email protected] (L.M.-S.) 2Network Research on Rare Diseases (CIBERER), U714, C/Melchor Fernández Almagro, 3, 28029 Madrid, Spain; [email protected] (S.L.); marta.carreter[email protected] (M.C.) 3Regenerative Medicine and Tissue Engineering Group, Health Research Institute-Jiménez Díaz Foundation University Hospital (IIS-FJD), Av. de los Reyes Católicos, 2, 28040 Madrid, Spain 4Centre for Energy, Environment and Technology Research (CIEMAT), Av. Complutense 40, 28040 Madrid, Spain 5Bioinformatics and Biostatistics Unit, Prince Felipe Research Center (CIPF), C/Eduardo Primo Yúfera, 3, 46012 Valencia, Spain; [email protected] 6Spanish National Bioinformatics Institute, ELIXIR-Spain (INB, ELIXIR-ES), C/Eduardo Primo Yúfera, 3, 46012 Valencia, Spain 7 Tissue Engineering Unit, Blood and Tissue Community Center of Asturias (CCST), C/Emilio Rodríguez Vigil s/n, 33006 Oviedo, Spain; [email protected] 8Clinical Bioinformatics Research Area, Progress and Health Foundation (FPS), CDCA, Hospital Virgen del Rocío, Av. Manuel Siurot s/n, 41013 Sevilla, Spain; [email protected] 9Functional Genomics Node, INB-ELIXIR-es, Progress and Health Foundation (FPS), Hospital Virgen del Rocío, Av. Manuel Siurot s/n, 41013 Sevilla, Spain 10 Institute of Biomedicine of Seville (IBIS), Hospital Virgen del Rocio, Av. Manuel Siurot s/n, 41013 Sevilla, Spain 11 Servicio de Análisis de Sistemas de Información Sanitaria, Consellería de Sanitat Universal i Salut Pública, Generalitat Valenciana, Carrer de Misser Mascó, 31-33, 46010 Valencia, Spain *Correspondence: [email protected]; Tel.: +34-91-624207 Abstract: Defective healing leading to cutaneous ulcer formation is one of the most feared complications of diabetes due to its consequences on patients’ quality of life and on the healthcare system. A more in-depth analysis of the underlying molecular pathophysiology is required to develop effective healing-promoting therapies for those patients. Major architectural and functional differences with human epidermis limit extrapolation of results coming from rodents and other small mammal-healing models. Therefore, the search for reliable humanized models has become mandatory. Previously, we developed a diabetes-induced delayed humanized wound healing model that faithfully recapitulated the major histological features of such skin repair-deficient condition. Herein, we present the results of a transcriptomic and functional enrichment analysis followed by a mechanistic analysis performed in such humanized wound healing model. The deregulation of genes implicated in functions such as angiogenesis, apoptosis, and inflammatory signaling processes were evidenced, confirming published data in diabetic patients that in fact might also underlie some of the histological features previously reported in the delayed skin-humanized healing model. Altogether, these molecular findings support the utility of such preclinical model as a valuable tool to gain insight into the molecular basis of the delayed diabetic healing with potential impact in the translational medicine field. Keywords: transcriptomics; wound healing; diabetes; skin-humanized mice; enrichment analysis Genes 2021,12, 47. https://doi.org/10.3390/genes12010047 https://www.mdpi.com/journal/genes
Genes 2021,12, 47 2 of 13 1. Introduction Diabetes is a systemic chronic disorder with a high and continuously increasing incidence and its global prevalence is projected to reach 10.2% by 2030 [ 1 ]. Impaired wound healing is one of the major complications associated to diabetes; in fact, between 19% and 34% of diabetic patients will develop a foot ulcer during the course of their illness [ 2 ]. High ulcer recurrence rates as well as frequent infections contribute to prolonged hospitalization and lead to an increased risk of lower limb amputation which has significant effects on patient morbidity and mortality [ 2 ]. Consequently, the quality of life of patients is often compromised. Furthermore, high costs associated with the clinical management of recurrent diabetic wounds represent a significant economic impact on the healthcare system. A persistent inflammatory response, impaired fibroblast function, hyperproliferative non-migratory epidermis, decreased angiogenesis, altered extracellular matrix (ECM) deposition, increased levels of proteases, and unresponsiveness to growth factor signals are the most recognized features of non-healing diabetic wounds [ 3 – 6 ]. However, the incomplete understanding of the underlying mechanisms responsible for defective diabetic healing contributes to the lack of effective treatments [ 7 ]. The study of diabetic wound healing in patients is limited mainly by technical and ethical considerations as well as by the complexity and heterogeneity of the disease. For these reasons, a plethora of murine models of diabetes have been developed [ 8 – 10 ] with the consequent flaws in extrapolation of the results due to functional and structural differences between human and rodent models. On the other hand, studies in large animals such as pigs are troublesome and expensive. In this context, the use of humanized mouse models poses a great potential for mimicking human conditions more accurately while preserving some of the advantages of the use of animal experimentation. Our group has extensive experience in generating humanized models based on the permanent engraftment of human bioengineered skin onto the back of immunodeficient mice [ 11 – 14 ] that in fact resemble a wide variety of physiological and pathological cutaneous processes, including wound healing [ 15 , 16 ]. Specifically, we developed a diabetesinduced delayed humanized wound healing model that reproduced some of the main cutaneous features described in the healing impairment of diabetic patients, including a prolonged inflammatory response, poor dermal matrix remodeling, impaired angiogenic response, and impaired migratory activity of keratinocytes [ 17 ]. These humanized models have also been suitable platforms to test pharmacologic, cell and gene therapeutic approaches [13,15–18]. A broad variety of analytical tools including standard cell, molecular, and different “omics” technologies have been used for unravelling deregulated processes in wound healing [ 19 – 21 ]. In the present work, a transcriptomic study together with a functional enrichment analysis of the wound healing response, followed by a more sophisticated mechanistic signaling pathway analysis [ 22 , 23 ], were performed in the diabetic skin-humanized mouse model aiming to compare our results with previous studies in diabetic patients and in such preclinical model. 2. Materials and Methods 2.1. Wound Healing Experimental Design in the Diabetic Skin Humanized Mouse Model NMRI nude mice (Rj:NMRI-Foxn1nu; Elevage Janvier Laboratories, Le Genest Saint Isle, France) were orthotopically grafted with bioengineered cutaneous equivalents [ 11 , 12 ]. This skin equivalent is based on human keratinocytes (epidermal component) seeded onto the top of the fibrin matrix populated with live human fibroblasts (dermal component). Cells from a skin biopsy of a healthy donor were isolated as previously described [ 11 , 12 , 24 ] after informed written consent and in accordance with the Helsinki declaration of 1975 and further revisions and with the Spanish regulation. The grafting of bioengineered human skin equivalents was performed under sterile conditions at the Centro de Investigaciones Energéticas Medioambientales y Tecnológicas (CIEMAT) Laboratory Animals
Genes 2021,12, 47 3 of 13 Facility (European registration number ES280790000183). All experimental procedures were performed in accordance with the corresponding regulations regarding experimental animal welfare. Experimental diabetes was induced in skin-humanized mice 10 weeks postgrafting by intraperitoneal injections of streptozotocin (STZ; Sigma-Aldrich, St. Luis, MO, USA), as previously described [ 17 ]. Wound healing experiments were performed in the stable human skin engrafted on immunodeficient mice exposed to sustained hyperglycemia for 6 weeks. Specifically, 2 mm-circular excisional wounds were created, and the excised tissue was harvested and used as a reference for gene expression analysis. After 24 h, a 6 mm-circular ring around the wound edge was taken. Samples were transferred immediately into liquid nitrogen for fast freezing and stored at − 80 ◦ C until processing. A total of 26 skin samples were obtained (n= 7 for control mice and n= 6 for diabetic mice for each condition, i.e., t = 0 h and t = 24 h). 2.2. RNA Extraction Total RNA isolation was performed using TRIzol ® (Invitrogen, Carlsbad, CA, USA) extraction-based methods, followed by purification in columns (Qiagen, Hilden, Germany). Briefly, skin samples were completely lysed in 1 mL of TRIzol, and subsequently 0.2 mL of chloroform was added to the suspension. After incubation at room temperature for 5 min to allow phase separation, the mixture was centrifuged at 14,000 rpm (15 min at 4 ◦ C) and RNA was isolated from the aqueous phase. RNA was then precipitated with isopropanol, washed with 70% ethanol, and resuspended in RNase-free water (Qiagen). Total isolated RNA was further purified with the RNeasy kit (Qiagen). RNA concentration was determined using a NanoDrop ™ . Spectrophotometer (Thermo Scientific, Waltham, MA, USA) and RNA integrity was verified with a Bioanalyzer (Agilent, Santa Clara, CA, USA). 2.3. Microarray Analysis Total RNA extracted from 26 skin samples was used to generate double-stranded cDNA. Briefly, biotin-labeled cRNA was fragmented and hybridized overnight to an Affymetrix Human Genome GeneChip ® (HU133 2.0 GeneChip, 22277 probesets), according to the manufacturer’s protocol. After 16 h of hybridization at 45 ◦ C, arrays were washed, stained with streptavidin-phycoerythrin, and then scanned with the GeneArray Scanner (Affymetrix, Santa Clara, CA, USA). The quality control of the scanned images included the examination of visible artifacts, the confirmation of proper grid alignment and the subtraction of background intensity. GeneChip ® operating software (GCOS) was used to generate “.CEL” files. All raw data comply with the Minimum Information About a Microarray Experiment (MIAME) guidelines. Gene expression DataSets are available in the Gene Expression Omnibus (GEO) on the NCBI website (http://www.ncbi.nlm.nih.gov/geo; accession number GSE147890). 2.4. Data Processing and Statistical Analysis Raw data from microarrays was normalized using the rma algorithm from the affy package (Bioconductor 3.12 Released) [ 25 ]. This algorithm includes background correction, quantile normalization and probeset summarization using the median. A Principal Component Analysis (PCA) was performed to the normalized data in order to detect the presence of outliers. Differential expression analysis between groups was evaluated with the limma package from Bioconductor [ 26 ] using the Babelomics suite v5.0 [ 27 ]. The p-values were corrected using the Benjamini–Hochberg method for multiple testing and FDR [ 28 ]. Finally, a functional enrichment analysis was performed using the updated versions of Gene Ontology [ 29 ] and the Kyoto Encyclopedia of Genes and Genomes (KEGG) [ 30 ] annotation databases through the tool incorporated in the Database for Annotation, Visualization, and Integrated Discovery (DAVID v6.8) [31].
Genes 2021,12, 47 4 of 13 2.5. Pathway Activity Analysis A signaling circuit activity analysis method implemented in the Hipathia tool [ 23 ] was employed to study wound healing in both control and diabetic mice. Under this approach, signaling circuits are defined within KEGG signaling pathways as the chain of proteins that connect a receptor protein to effector proteins that trigger specific cellular activities. Normalized gene expression values are taken as proxies of protein activity. 3. Results 3.1. Exploratory Microarray Data Analysis RNA extracts of 26 skin-humanized samples were hybridized to Affymetrix GeneChip ® Human Genome U133 2.0 Array. The total expression of >20,000 probes was normalized and preprocessed. The normalized data for all probes in all samples is provided in Table S1 . First, in order to ensure the highest data quality, a PCA was performed. Two control samples taken at 0 h were found as outliers (95% confidence) and discarded to avoid technical variability bias in the data analysis. The remaining 24 samples were clearly grouped into 4 different categories: control samples at 0 h and 24 h, and diabetic samples at 0 h and 24 h (Figure 1). Specifically, the bigger difference in gene expression between groups was related to wound healing process (0 h versus 24 h). In fact, PC1 mainly separated samples from 0 h and 24 h (57.3% variance) while PC2 (6.8% variance) was the main factor discriminating between the diabetic and control mice (Figure 1). Genes2021,12,xFORPEERREVIEW4of14 [30]annotationdatabasesthroughthetoolincorporatedintheDatabaseforAnnotation, Visualization,andIntegratedDiscovery(DAVIDv6.8)[31]. 2.5.PathwayActivityAnalysis AsignalingcircuitactivityanalysismethodimplementedintheHipathiatool[23] wasemployedtostudywoundhealinginbothcontrolanddiabeticmice.Underthisap‐ proach,signalingcircuitsaredefinedwithinKEGGsignalingpathwaysasthechainof proteinsthatconnectareceptorproteintoeffectorproteinsthattriggerspecificcellular activities.Normalizedgeneexpressionvaluesaretakenasproxiesofproteinactivity. 3.Results 3.1.ExploratoryMicroarrayDataAnalysis RNAextractsof26skin‐humanizedsampleswerehybridizedtoAffymetrix GeneChip ® HumanGenomeU1332.0Array.Thetotalexpressionof>20,000probeswas normalizedandpreprocessed.Thenormalizeddataforallprobesinallsamplesispro‐ videdinTableS1.First,inordertoensurethehighestdataquality,aPCAwasperformed. Twocontrolsamplestakenat0hwerefoundasoutliers(95%confidence)anddiscarded toavoidtechnicalvariabilitybiasinthedataanalysis.Theremaining24sampleswere clearlygroupedinto4differentcategories:controlsamplesat0hand24h,anddiabetic samplesat0hand24h(Figure1).Specifically,thebiggerdifferenceingeneexpression betweengroupswasrelatedtowoundhealingprocess(0hversus24h).Infact,PC1 mainlyseparatedsamplesfrom0hand24h(57.3%variance)whilePC2(6.8%variance) wasthemainfactordiscriminatingbetweenthediabeticandcontrolmice(Figure1). Figure1.PrincipalComponentAnalysis(PCA)biplotofallnormalizedsamples(n=24)fromthe4 groups(C0:controlsamplesat0h;C24:controlsamplesat24h;D0:diabeticsamplesat0h;D24: diabeticsamplesat24h).Circlesshow95%confidence. DifferentialexpressionanalysiswasperformedusinglimmapackagefromBiocon‐ ductor.AsinthePCAanalysis,thenumberofdifferentiallyexpressed(DE)probesindi‐ catedthatthemaindrivingforceintermsofdifferentialgeneexpressionwasmainlytrig‐ geredinresponsetothecutaneouswoundduringthefirst24h(Table1).Indeed,and Figure 1. Principal Component Analysis (PCA) biplot of all normalized samples (n= 24) from the 4 groups (C0: control samples at 0 h; C24: control samples at 24 h; D0: diabetic samples at 0 h; D24: diabetic samples at 24 h). Circles show 95% confidence. Differential expression analysis was performed using limma package from Bioconductor. As in the PCA analysis, the number of differentially expressed (DE) probes indicated that the main driving force in terms of differential gene expression was mainly triggered in response to the cutaneous wound during the first 24 h (Table 1). Indeed, and since the experimental diabetes process was sustained for a longer period (6 weeks), a minor number of genes was affected at the time of sampling.
Genes 2021,12, 47 5 of 13 Table 1. Number of probes differentially regulated (FDR < 0.05) in each group comparison and their corresponding enriched Gene Ontology Biological Processes (GO_BPs) and KEGG pathways. All the differentially expressed probes in each comparison together with their p-value and fold change can be found in Table S2. Comparison Differentially Expressed Probes Up-Regulated Probes Down-Regulated Probes Significant GO_BPs Significant KEGG C24 vs. C0 7570 3350 4220 456 42 D0 vs. C0 403 101 302 100 7 D24 vs. C24 49 14 35 3 0 D24 vs. D0 8686 3801 4885 519 71 3.2. Differential Expression and Functional Enrichment Analysis of the Diabetes-Induced Skin-Humanized Mouse Model In relation to the experimentally induced diabetes in the skin-humanized mouse model, 403 probes were differentially expressed (FDR < 0.05) between both groups of animals before wounding (D0 vs. C0), mainly in response to the diabetes induction. Likewise, the number of biological processes (GO_BPs) and KEGG pathways associated with those dysregulated probes followed the same trend (Table 1). However, the number of DE probes between diabetic and control mice at 24 h (D24 vs. C24) was much lower (49 probes). This could be explained due to the bigger variability of samples taken at 24 h compared to the samples at 0 h which in fact, formed a more compact cluster (Figure 1). Enrichment analysis of those 403 DE probes using DAVID bioinformatics tool showed a total of 100 biological processes from the Gene Ontology and 7 KEGG pathways ( FDR < 0.05 ) associated to the gene expression changes during the process of experimental diabetes induction (i.e., D0 vs C0 comparison) (Table S3). Specifically, terms such as ECM organization, cell-cell adhesion as well as PI3K/Akt signaling, ECM receptor interaction and focal adhesion were altered in the diabetes induction process (Figure 2). Genes2021,12,xFORPEERREVIEW5of14 sincetheexperimentaldiabetesprocesswassustainedforalongerperiod(6weeks),a minornumberofgeneswasaffectedatthetimeofsampling. Table1.Numberofprobesdifferentiallyregulated(FDR<0.05)ineachgroupcomparisonand theircorrespondingenrichedGeneOntologyBiologicalProcesses(GO_BPs)andKEGGpathways. Allthedifferentiallyexpressedprobesineachcomparisontogetherwiththeirp‐valueandfold changecanbefoundinTableS2. Compari‐ son DifferentiallyEx‐ pressedProbes Up‐Regu‐ latedProbes Down‐Regu‐ latedProbes Significant GO_BPs Significant KEGG C24vs.C075703350422045642 D0vs.C04031013021007 D24vs.C2449143530 D24vs.D086863801488551971 3.2.DifferentialExpressionandFunctionalEnrichmentAnalysisoftheDiabetes‐InducedSkin‐ HumanizedMouseModel Inrelationtotheexperimentallyinduceddiabetesintheskin‐humanizedmouse model,403probesweredifferentiallyexpressed(FDR<0.05)betweenbothgroupsofan‐ imalsbeforewounding(D0vs.C0),mainlyinresponsetothediabetesinduction.Like‐ wise,thenumberofbiologicalprocesses(GO_BPs)andKEGGpathwaysassociatedwith thosedysregulatedprobesfollowedthesametrend(Table1).However,thenumberofDE probesbetweendiabeticandcontrolmiceat24h(D24vs.C24)wasmuchlower(49 probes).Thiscouldbeexplainedduetothebiggervariabilityofsamplestakenat24h comparedtothesamplesat0hwhichinfact,formedamorecompactcluster(Figure1). Enrichmentanalysisofthose403DEprobesusingDAVIDbioinformaticstool showedatotalof100biologicalprocessesfromtheGeneOntologyand7KEGGpathways (FDR<0.05)associatedtothegeneexpressionchangesduringtheprocessofexperimental diabetesinduction(i.e.,D0vsC0comparison)(TableS3).Specifically,termssuchasECM organization,cell‐celladhesionaswellasPI3K/Aktsignaling,ECMreceptorinteraction andfocaladhesionwerealteredinthediabetesinductionprocess(Figure2). Figure2.BargraphshowingthedysregulatedGeneOntologybiologicalprocesses(GO_BPs)andKEGGpathways(FDR <0.05)withthelowestp‐value(p˂ 0.01)inthediabetesinductionprocess(D0vsC0comparison).Yellowbarsindicatethe p‐valueoftheterminlogformandorangebarsindicatethepercentageofdysregulatedgenesineachterm. 3.3.DifferentialExpressionandFunctionalEnrichmentAnalysisoftheWoundHealingProcess intheDiabetes‐InducedSkin‐HumanizedMouseModel Figure 2. Bar graph showing the dysregulated Gene Ontology biological processes (GO_BPs) and KEGG pathways ( FDR < 0.05 ) with the lowest p-value (p< 0.01) in the diabetes induction process (D0 vs C0 comparison). Yellow bars indicate the p-value of the term in log form and orange bars indicate the percentage of dysregulated genes in each term. 3.3. Differential Expression and Functional Enrichment Analysis of the Wound Healing Process in the Diabetes-Induced Skin-Humanized Mouse Model Aiming to study the functional similarities and differences in both wound healing processes (C24vsC0 and D24vsD0), two enrichment analyses were carried out. Specifically, the common dysregulated probes and also the specific probes to each process were identified (Figure 3).
Genes 2021,12, 47 6 of 13 Genes2021,12,xFORPEERREVIEW6of14 Aimingtostudythefunctionalsimilaritiesanddifferencesinbothwoundhealing processes(C24vsC0andD24vsD0),twoenrichmentanalyseswerecarriedout.Specifi‐ cally,thecommondysregulatedprobesandalsothespecificprobestoeachprocesswere identified(Figure3). Figure3.Venndiagramofcommonanduniquedifferentiallyexpressedprobesinbothwoundhealingprocesses(C24vsC0 andD24vsD0).Thenumberineachbargraphindicatestheup‐(red)anddown‐regulated(blue)probes(FDR<0.05).(*): Sevenprobesetsoutofthetotal5902commononesshowedoppositeregulation. 3.3.1.CommonTranscriptomicResponsetoWoundHealinginBothControlandDia‐ beticMice Outofthe7570and8686probesdysregulatedinthewoundhealingprocessofcon‐ trolanddiabeticmicerespectively(C24vsC0andD24vsD0;Figure3),5902probeswere commoninbothgroups,andthereforelinkedtothesamedysregulatedfunctionsand pathways.Thefoldchangeoftheseprobeswasalsosurprisinglysimilarinbothcompar‐ isons(FigureS1A),displayingastronglinearcorrelation(slope=1.074;R2=0.93).How‐ ever,sevenprobesoutofthose5902hadadifferentexpressionpatternbetweengroups (i.e.,upregulatedincontrolsanddownregulatedindiabeticanimals,orviceversa) (FigureS1B)whichsuggestsanoppositemechanisminthefunctionsrelatedtothose sevengenesinthedifferentexperimentalgroups. Anenrichmentanalysisofthecommontranscriptomicresponsetowounding(5895 probesdisplayingthesameexpressionpattern)inbothcontrolanddiabeticmicewascar‐ riedout.Thus,upregulatedanddownregulatedprobesat24hpost‐wounding(2435and 3460probesrespectively;Figure3)wereanalyzedseparately.Specifically,downregulated geneswerepredominantlyenrichedinfunctionsrelatedtotranscriptionprocesses,gene regulationandECMorganization(Figure4A,B).Ontheotherhand,upregulatedprobes wereenrichedinaverywidearrayoffunctions,manyofthemrelatedtotheinflammatory responseafterwounding,suchasWntsignalingpathway,NFkβsignaling,TNFmediated Figure 3. Venn diagram of common and unique differentially expressed probes in both wound healing processes (C24vsC0 and D24vsD0). The number in each bar graph indicates the up-(red) and down-regulated (blue) probes (FDR < 0.05). (*): Seven probesets out of the total 5902 common ones showed opposite regulation. 3.3.1. Common Transcriptomic Response to Wound Healing in Both Control and Diabetic Mice Out of the 7570 and 8686 probes dysregulated in the wound healing process of control and diabetic mice respectively (C24vsC0 and D24vsD0; Figure 3), 5902 probes were common in both groups, and therefore linked to the same dysregulated functions and pathways. The fold change of these probes was also surprisingly similar in both comparisons (Figure S1A), displaying a strong linear correlation (slope = 1.074; R2 = 0.93). However, seven probes out of those 5902 had a different expression pattern between groups (i.e., upregulated in controls and downregulated in diabetic animals, or vice versa) (Figure S1B) which suggests an opposite mechanism in the functions related to those seven genes in the different experimental groups. An enrichment analysis of the common transcriptomic response to wounding (5895 probes displaying the same expression pattern) in both control and diabetic mice was carried out. Thus, upregulated and downregulated probes at 24 h post-wounding (2435 and 3460 probes respectively; Figure 3) were analyzed separately. Specifically, downregulated genes were predominantly enriched in functions related to transcription processes, gene regulation and ECM organization (Figure 4A,B). On the other hand, upregulated probes were enriched in a very wide array of functions, many of them related to the inflammatory response after wounding, such as Wnt signaling pathway, NFkβsignaling, TNF mediated response, as well as DNA replication, cell–cell adhesion or MAPK signaling pathways and angiogenesis, among others (Figure 4C,D).
Genes 2021,12, 47 7 of 13 Genes2021,12,xFORPEERREVIEW7of14 response,aswellasDNAreplication,cell–celladhesionorMAPKsignalingpathwaysand angiogenesis,amongothers(Figure4C,D). Figure4.Diagramofcommonenrichedfunctions(GO_BPs:GeneOntologybiologicalprocesses,andKEGGpathways)in thedownregulated(A,B)andupregulatedprobes(C,D)inresponsetowoundinginbothgroups(C24vsC0andD24vsD0). ThecompletelistofenrichedGO_BPsandKEGGpathwaysisavailableinTableS4(sheettabsAandB). 3.3.2.DifferentTranscriptomicResponsetoWoundHealinginControlandDiabetic Mice Afunctionalenrichmentanalysisoftheprobeswhichalteredexpressionwasunique toeithercontrol(1668)ordiabetichealing(2784)(Figure3),wasperformed(Table2).Spe‐ cifically,upregulatedanddownregulatedprobeswerestudiedseparately,aimingtoun‐ derstandthewoundhealingprocessinherenttoeachexperimentalgroup. Table2.NumberofsignificantGeneOntologybiologicalprocesses(GO_BPs)andKEGGpath‐ waysenrichedintheuniqueup‐anddown‐regulatedprobes(FDR<0.05)inbothC24vsC0and D24vsD0comparisons.ThelistofthesignificanttermscanbefoundinTableS4(sheettabsC–F). Comparison SignificantGO_BPsSignificantKEGG Up‐Regulated Probes Down‐Regu‐ latedProbes Up‐Regulated Probes Down‐Regu‐ latedProbes UniqueC24vsC0671221127 UniqueD24vsD020486728 Interestingly,amarkedenrichmentforGO_BPsandKEGGpathways(204and72 functions,respectively)wasobservedinupregulatedprobesinthediabeticwoundheal‐ ingprocess(uniqueinD24vsD0).Specifically,pathwaysrelatedtocellularmovementand responsetomechanicalstimulusaswellastotheinflammatoryresponse,suchasTGFβ, NFkβ,TNF,andchemokinesignalingpathways,wereenrichedinthoseprobes (TableS4C).TherewerealsomanyuniqueGO_BPsrelatedtoderegulationofapoptosis processes. Figure 4. Diagram of common enriched functions (GO_BPs: Gene Ontology biological processes, and KEGG pathways) in the downregulated ( A , B ) and upregulated probes ( C , D ) in response to wounding in both groups (C24vsC0 and D24vsD0). The complete list of enriched GO_BPs and KEGG pathways is available in Table S4 (sheet tabs A and B). 3.3.2. Different Transcriptomic Response to Wound Healing in Control and Diabetic Mice A functional enrichment analysis of the probes which altered expression was unique to either control (1668) or diabetic healing (2784) (Figure 3), was performed (Table 2). Specifically, upregulated and downregulated probes were studied separately, aiming to understand the wound healing process inherent to each experimental group. Table 2. Number of significant Gene Ontology biological processes (GO_BPs) and KEGG pathways enriched in the unique upand down-regulated probes (FDR < 0.05) in both C24vsC0 and D24vsD0 comparisons. The list of the significant terms can be found in Table S4 (sheet tabs C–F). Comparison Significant GO_BPs Significant KEGG Up-Regulated Probes Down-Regulated Probes Up-Regulated Probes Down-Regulated Probes Unique C24vsC0 67 122 11 27 Unique D24vsD0 204 86 72 8 Interestingly, a marked enrichment for GO_BPs and KEGG pathways (204 and 72 functions, respectively) was observed in upregulated probes in the diabetic wound healing process (unique in D24vsD0). Specifically, pathways related to cellular movement and response to mechanical stimulus as well as to the inflammatory response, such as TGF β , NFk β , TNF, and chemokine signaling pathways, were enriched in those probes (Table S4C). There were also many unique GO_BPs related to deregulation of apoptosis processes. On the other hand, unique pathways in the control wound healing process (C24vsC0) were mainly related to mitochondrial and oxidative phosphorylation processes (upregulated genes), as well as to transcriptional regulation (downregulated genes) (Table S4, sheet tabs E and F).
Genes 2021,12, 47 8 of 13 Aiming to identify processes equally altered in both wound healing processes but due to the dysregulation of different genes, coincident and discordant functions enriched in unique probes in C24vsC0 and D24vsD0 were evaluated (Figure 5). Genes2021,12,xFORPEERREVIEW8of14 Ontheotherhand,uniquepathwaysinthecontrolwoundhealingprocess(C24vsC0) weremainlyrelatedtomitochondrialandoxidativephosphorylationprocesses(upregu‐ latedgenes),aswellastotranscriptionalregulation(downregulatedgenes)(TableS4, sheettabsEandF). Aimingtoidentifyprocessesequallyalteredinbothwoundhealingprocessesbut duetothedysregulationofdifferentgenes,coincidentanddiscordantfunctionsenriched inuniqueprobesinC24vsC0andD24vsD0wereevaluated(Figure5). Figure5.Venndiagramofthenumberoffunctions(GeneOntologybiologicalprocessesand KEGGpathways)alteredinuniqueprobesupanddownregulatedinC24vsC0(1668DEprobes) andD24vsD0(2784DEprobes)groups,respectively.Whitecircleshowscoincidentfunctionsal‐ teredbygeneswithoppositeregulation(upregulatedinD24vsD0anddownregulatedin C24vsC0). Itmustbenotedthatthenumberofcoincidentpathwayswassignificantlyhigher betweenupregulatedgenesinthediabeticwoundhealingandthedownregulatedones incontrolsamples(asindicatedbyawhitecircleinFigure5andhighlightedinyellowin TableS4).Thisfactsuggeststhatseveralfunctionsormetabolicrouteswerebeingtrig‐ geredinoppositedirectionsinresponsetowoundhealinginbothexperimentalgroups, includingregulationofcellproliferation,fibroblastgrowthfactorsignalingpathways, FoxOsignalingpathway,TFGβreceptorsignalingpathway,cellmigration,andECMor‐ ganization,amongothers(TableS4).Someofthesefunctionshavebeenalreadyassociated toimpaireddiabetichealing[5,32,33]. 3.4.MechanisticSignalingPathwayAnalysisoftheWoundHealing Tofurtherdescribetheunderlyingmolecularmechanismsinvolvedinthewound healingprocess,amechanisticactivityanalysisthatdecomposesKEGGsignalingpath‐ waysintocanonicalsub‐pathways(circuits)wasemployed.Thismethodconnectsrecep‐ torproteins,bymeansofachainofintermediateproteins,toeffectorproteinsthattrigger specificcellactivities.Thisanalysiswasperformedusingthenormalizedexpressionval‐ uesofbothC24vsC0andD24vsD0comparisons.Allsignificantcircuits(FDR<0.05)and effectorproteinsareincludedinTableS5.Specifically,atotalof382circuitsinvolvedin Figure 5. Venn diagram of the number of functions (Gene Ontology biological processes and KEGG pathways) altered in unique probes up and downregulated in C24vsC0 (1668 DE probes) and D24vsD0 (2784 DE probes) groups, respectively. White circle shows coincident functions altered by genes with opposite regulation (upregulated in D24vsD0 and downregulated in C24vsC0). It must be noted that the number of coincident pathways was significantly higher between upregulated genes in the diabetic wound healing and the downregulated ones in control samples (as indicated by a white circle in Figure 5and highlighted in yellow in Table S4). This fact suggests that several functions or metabolic routes were being triggered in opposite directions in response to wound healing in both experimental groups, including regulation of cell proliferation, fibroblast growth factor signaling pathways, FoxO signaling pathway, TFG β receptor signaling pathway, cell migration, and ECM organization, among others (Table S4). Some of these functions have been already associated to impaired diabetic healing [5,32,33]. 3.4. Mechanistic Signaling Pathway Analysis of the Wound Healing To further describe the underlying molecular mechanisms involved in the wound healing process, a mechanistic activity analysis that decomposes KEGG signaling pathways into canonical sub-pathways (circuits) was employed. This method connects receptor proteins, by means of a chain of intermediate proteins, to effector proteins that trigger specific cell activities. This analysis was performed using the normalized expression values of both C24vsC0 and D24vsD0 comparisons. All significant circuits (FDR < 0.05) and effector proteins are included in Table S5. Specifically, a total of 382 circuits involved in 53 KEGG pathways (D24vsD0) and 424 circuits in 55 pathways (C24vsC0) were affected in each comparison. The analysis of the common and unique altered circuits evidenced that most of them (262 circuits) were common to both wound healing processes while only 102 and 137 circuits were unique in D24vsD0 and C24vsC0, respectively. In both cases, these circuits were mainly involved in KEGG pathways already highlighted in the enrichment analysis such as TGF β , TNF, NFk β , PI3K, Wnt, and FoxO signaling pathways, as well
Genes 2021,12, 47 9 of 13 as apoptotic processes. In addition, there was a total coincidence in the direction (up or downregulation) of the common circuits between both comparisons. 4. Discussion There is still a lack of understanding of the underlying diabetes wound repair impairment. Since the use of human samples has ethical and technical issues that cannot be easily ignored [ 21 ], different animal models of diabetes have been developed in the last decades [ 8 – 10 ]. Although a huge amount of valuable basic knowledge has been obtained from such small pre-clinical animal models, they do not totally reproduce the complexity of the human disease, thus it is mandatory to recreate humanized contexts aiming to shed light on the molecular and cellular mechanisms responsible for diabetes–associated impaired healing. The humanized mouse model developed by our group [ 11 , 12 , 24 ] has been proven as a suitable platform to perform wound-healing studies in a humanized context since all major features of human cutaneous wound healing are accurately recapitulated [ 15 – 17 ]. In addition, we have developed humanized models that recreate a wide range of human cutaneous disorders in which interactions between the immune system, epidermis and the environment likely occur such as psoriasis and different genodermatosis [ 13 ]. Currently, humanized mice carrying both a functioning immune system and skin of human origin are being developed aimed to better recapitulate human cutaneous physiology [ 34 ]. These multi-tissue compartment models, though promising, need further research. Herein, we present the transcriptomic analysis of a humanized model of delayed wound healing in a diabetic context in order to dissect the early molecular mechanisms underlying diabetes wound repair impairment in such preclinical model. Since the major transcriptomic changes are usually observed within the first 24 h after wounding [ 35 , 36 ], the samples were taken at this time point. An enrichment and a mechanistic pathway analysis were also performed using the Gene Ontology and KEGG databases, thus the comparison was done at a functional level, and not only based in lists of dysregulated genes. Furthermore, this approach also allowed us to compare biological processes or pathways that may be equally altered in both diabetic and control wound healing processes, although through different gene expression variations. The analysis of the transcriptomic data in the chemically induced diabetic skinhumanized mouse model showed that there were more differences in the gene expression due to the diabetes induction at the basal level (D0 vs. C0) than 24 h after wounding (D24 vs. C24). Indeed, this observation might be due to an acute genetic response to healing at 24 h in both experimental groups that would be drastically modifying the global gene expression and thus, minimizing the intrinsic differences between diabetic and control mice. This fact could also explain that a minor number of probes were deregulated in unwounded skin samples (D0 vs. C0) than in both healing processes at 24 h (D24 vs. D0 and C24 vs. C0). In general terms, the differential global gene expression in the skin of the diabetic humanized mice (D0vsC0) closely reproduced the cutaneous transcriptome profiling in diabetic patients [ 37 , 38 ]. Specifically, those studies showed altered GO_BPs and KEGG pathways related to focal adhesions, cell adhesions mediated by integrins, and ECMrelated processes, as in our diabetic skin-humanized mouse model (Figure 2). In addition, other relevant functions such as PI3K/Akt signaling pathway and collagen deposition mechanisms reported in diabetic patients [ 39 , 40 ] were also identified in our humanized model. Nevertheless, it must be noted that, unlike long-term exposure to hyperglycemia in diabetic patients, samples in the diabetic skin-humanized mouse model were taken after six weeks of sustained high glucose levels. Therefore, other deregulated molecular mechanisms triggered later by longer exposition to hyperglycemia might also be involved. Notwithstanding this fact, the main cutaneous features reported in diabetic patients, such as decreased innervation and vascularization, were faithfully mimicked in the regenerated human skin of diabetic mice after six weeks of sustained high glucose levels [17].