Abstract

   Mechanical ventilation contributes to diaphragm atrophy and muscle
   weakness, which is referred to as ventilator-induced diaphragmatic
   dysfunction (VIDD). The pathogenesis of VIDD has not been fully
   understood until recently. The aim of this study was to investigate the
   effects of 24 h of mechanical ventilation on fibro-adipogenic
   progenitor (FAP) proliferation, endothelial-mesenchymal transition
   (EndMT), and immune cell infiltration driving diaphragm fibrosis in a
   rabbit model. The rabbits were anaesthetized and randomly divided into
   two groups (n = 3 each group): a control group and an experimental
   group. Diaphragm nuclei for sequencing were prepared by dissociating
   and filtering muscle tissue. 10X Genomics Platform for single-nucleus
   RNA sequencing (snRNA-seq) was used to profile the cells. Normalization
   and clustering were performed by Seurat, and clusters were manually
   annotated as different cell types. In this study, we performed
   differentially expressed genes (DEGs) analysis, Kyoto Encyclopedia of
   Genes and Genomes (KEGG) analysis, pseudotime analysis and high
   dimensional weighted gene coexpression network analysis (hdWGCNA) to
   identify the key genes and signaling pathways related to the
   pathogenesis of VIDD. We further performed quantitative real-time
   polymerase chain reaction (qRT-PCR) and Western blotting to verify the
   results of snRNA-seq. The snRNA-seq results showed that acute
   postmechanical ventilation diaphragm cell changes included an increase
   in the proportion of fibroblasts and a decrease in the proportion of
   myofibres. The DEGs, KEGG, hdWGCNA and pseudotime analyses demonstrated
   that fibro-adipogenic progenitor (FAP) proliferation,
   endothelial-mesenchymal transition (EndMT) and immune cell infiltration
   are the three main processes involved in early stage of fibrosis
   development, among which Pdgfd, Sema3a, Cxcr2, are the corresponding
   regulatory genes. Glycolysis and the gene Pfkfb3 are also important
   metabolic factors for fibrosis formation. Negr1 and Mef2c are involved
   in phrenic nerve ending loss and diaphragm fibre atrophy. The qRT-PCR
   data showed that the mRNA levels of the genes Pdgfd, Cxcr2, Pfkfb3 and
   Negr1 were significantly greater in the experimental group than in the
   control group (P < 0.01), and the expression levels of Sema3a and Mef2c
   were significantly lower (P < 0.01). Despite limitations, including the
   lack of functional evaluations to confirm ventilator-induced diaphragm
   dysfunction (VIDD) and the absence of data validating diaphragm
   unloading during ventilation, our findings suggest that FAP
   proliferation and immune cell infiltration may play a role in the early
   stage of driving diaphragm fibrosis during mechanical ventilation.
   However, future studies are needed to confirm these findings and
   investigate the potential mechanisms underlying them.

Supplementary Information

   The online version contains supplementary material available at
   10.1038/s41598-024-82530-4.

   Keywords: Single-nucleus RNA sequencing, Diaphragm fibrosis,
   Fibro-adipogenic progenitor, High dimensional weighted gene
   coexpression network analysis, Endothelial-mesenchymal transition

   Subject terms: Molecular medicine, Diagnostic markers

Introduction

   Invasive mechanical ventilation (MV) was commonly used to treat the
   increased number of patients with acute respiratory distress syndrome
   during the COVID-19 pandemic in the ICU^[34]1. MV induces diaphragmatic
   muscle fibre atrophy and profound muscle weakness, termed
   ventilator-induced diaphragmatic dysfunction (VIDD), which results in
   delayed weaning and substantial morbidity and mortality. Reactive
   oxygen species (ROS) activate proteolytic muscle proteases, including
   caspases and calpains, deregulate the ubiquitin proteasome system and
   enhance autophagy, which is linked to a decrease in diaphragm protein
   synthesis and accelerated proteolysis, which are major contributors to
   VIDD^[35]2. Bruells’s results indicated that diaphragm muscle fibre
   size increased and contractile function recovered following 12 h of
   spontaneous breathing as a treatment for VIDD. The recovery of the
   diaphragm allows the re-establishment of cellular homeostasis of
   protein synthesis and breakdown^[36]3.

   In the context of Duchenne muscular dystrophy, fibrosis emerges as a
   pathological repair process that occurs more frequently than normal
   diaphragm regeneration, and this fibrosis directly impacts muscle force
   and physical performance, contributing to the progressive decline in
   muscle function characteristic of the disease^[37]4. In ventilated
   critically ill patients with COVID-19, Shi observed for the first time
   that fibrosis quickly develops in the diaphragm compared with that in
   patients without COVID-19^[38]5. Qian demonstrated that positive
   end-expiratory pressure (PEEP) application for 48 h causes collagen
   deposition and fibrosis in the diaphragm of mechanically ventilated
   rabbits^[39]6. Shi observed that ECM partially replaces the area lost
   by myofiber atrophy to provide structural stability of the myofibers
   and fibrosis quickly develops in the diaphragm of ventilated critically
   ill patients during ICU stay^[40]7. These findings suggest that
   diaphragmatic fibrosis is a complex pathological process that can occur
   in various situations and there are the potential link between VIDD and
   fibrosis which may lead to diaphragmatic stiffness and decreased
   contractility.

   The diaphragm is a complex heterogeneous tissue composed of
   multinucleated myofibres and mononuclear cells, including immune cells,
   endothelial cells, muscle stem cells (satellite cells), and nonmyogenic
   mesenchymal progenitor (e.g., fibro-adipogenic progenitor, or
   FAP)^[41]8. However, little is known about the mechanism by which the
   diaphragm recovers from VIDD and the response of cell populations such
   as satellite cells and FAP to injury. High-throughput spatial
   transcriptomics of mouse diaphragm cells revealed the activation of
   Myogenin 1 (Myod1)-expressing satellite cells near the injury site 2
   days after injury^[42]9. FAP expressing the markers platelet-derived
   growth factor receptor alpha (Pdgfrα) and Cd34 are interstitial
   populations of mesenchymal cells that can differentiate into
   fibroblasts and adipocytes^[43]10. Fibroblasts are the principal cell
   types that regulate ECM formation and the process of fibrosis. In
   response to tissue injury caused by mechanical ventilation, activated
   FAP proliferate more quickly than satellite cells, which may contribute
   to diaphragm fibrosis^[44]11.

   Transcriptome studies on skeletal muscle were conducted using bulk
   RNA-seq and high-throughput single-cell or snRNA-seq. Bulk RNA-seq
   provides only averaged gene expression profiles, which inhibits the
   ability to identify distinct cell types and their dynamics. Single-cell
   RNA sequencing (scRNA-seq) has been widely utilized to explore
   different cell types and developmental trajectories based on the
   diversity of cellular states and to broaden our understanding of
   cellular heterogeneity in skeletal muscle^[45]12,[46]13. Skeletal
   muscle is unique in its cellular composition because it is composed
   primarily of multinucleated myofibres, which are not easily isolated
   via single-cell dissociation due to their large size and elongated
   shape. SnRNA-seq has been applied to multinucleated myofibres for
   transcriptome studies as an alternative to scRNA-seq^[47]14,[48]15.
   Weighted gene coexpression network analysis (WGCNA) is a biological
   data analysis tool for constructing gene coexpression networks and
   identifying core genes that are considered potential disease drivers
   and are suitable for bulk RNA sequencing data. High-dimensional
   weighted gene coexpression network analysis (hdWGCNA) is a tool that
   offers various functionalities for identifying robust modules and hub
   genes and deciphering their associated biological functions across
   multiple cell scales and was more suitable to analyse the snRNA-seq
   data in this study^[49]16,[50]17.

   The cellular and molecular mechanisms of diaphragm fibrosis and atrophy
   during mechanical ventilation are unclear. In this study, we profiled
   the transcriptomic landscape of the diaphragm via snRNA-seq on the 10X
   Genomics Platform and performed hdWGCNA after one day of mechanical
   ventilation to determine the major diaphragm cell types, gene
   regulatory networks and signaling pathways changes. We aimed to
   evaluate the effects of 24 h of mechanical ventilation on
   fibro-adipogenic progenitor (FAP) proliferation,
   endothelial-mesenchymal transition (EndMT), and immune cell
   infiltration driving diaphragm fibrosis in rabbits. By using snRNA-seq,
   we hoped to gain insights into the cellular and molecular mechanisms
   underlying the early changes that occur in the diaphragm during
   mechanical ventilation.

Methods

Animals

   The rabbits were purchased from Shanxi Medical University. All animals
   weighing between 2.0 and 2.5 kg had ad libitum access to food and
   water. The rabbits were housed at 22 °C with a humidity of 30–70% under
   a 12 h light/12 h dark cycle. All animal care and experimental
   procedures associated with this research were approved by the
   Institutional Animal Care and Use Committees at Changzhi Medical
   College, and all experiments were performed in accordance with ARRIVE
   guidelines.

Experimental protocols

   The rabbits were anaesthetized with pentobarbital sodium
   (intraperitoneal injection, 50 mg/kg). After induction of anaesthesia,
   the rabbits were randomly divided into two groups (n = 3 each group): a
   control group with only anaesthetized spontaneously breathing animals
   and an experimental group with anaesthetized and mechanically
   ventilated ones. The ventilator (Harvard Apparatus, Holliston, MA, USA)
   was used in pressure-controlled mode with a PEEP of 3 cm H[2]O, a
   respiratory frequency of 40 breaths/min and an inspiratory: expiratory
   ratio of 1:2. At 12 and 24 h, blood gas analysis was performed for all
   the animals. The parameters were adjusted to maintain the normal range
   of PaO[2] and PaCO[2]. An arterial catheter was inserted into the
   femoral artery to measure arterial blood pressure. Body temperature was
   maintained at 37 °C using a temperature-controlled heating plate. After
   24 h of mechanical ventilation, the animals were euthanized, and
   diaphragm samples were removed. The tissues were flash-frozen in liquid
   nitrogen and stored at -80 °C until sufficient nuclei were collected
   for sequencing. In previous studies we already successfully established
   a 24-hour mechanical ventilation rabbit diaphragm rest animal
   model^[51]18.

Single-nucleus RNA sequencing

   We used snRNA-seq to analyse rabbit diaphragm cells from the control
   (n = 2) and MV (n = 2) groups (OE Biotech Co., Ltd. Shanghai, China).
   Nuclei were isolated from four samples of dissociated diaphragm tissue
   (10–15 mg). The cells were placed in lysis buffer containing 0.1% NP40,
   10 mM Tris HCl, 146 mM NaCl, 1 mM CaCl[2], 21 mM MgCl[2], and 40 U/mL
   RNase inhibitor for 7 min. Complete nuclear lysis was confirmed using
   trypan blue staining, and 1 ml of ST wash buffer was added. The samples
   were then filtered through a 40 μm cell strainer and centrifuged (500×g
   for 10 min at 4 °C) to form a crude pellet. The resultant nuclear
   suspension was pelleted in 100 µl of PBS and 1.0% BSA with a light
   centrifugation step (250×g for 5 min at 4 °C). Nuclei were counted
   under a microscope. The cell nuclei suspension was diluted with
   PBS + 1% BSA to a concentration of 700–1200/µL to reach the optimal
   range for loading on the 10X Genomics Chromium Chip. The concentration
   at which nuclei are loaded onto the 10X Genomics Platform is an
   important parameter affecting the number of nuclei available for
   downstream analysis. The nuclei were then loaded into the 10X Genomics
   Chromium Chip using Next GEM Single Cell 3ʹ Kits v3.1 (1000268) to
   perform cDNA library amplification. cDNA libraries were constructed
   using a Chromium Next GEM Chip G Single Cell Kit (1000020), and cDNA
   quality was evaluated by fragment analysis (5200 Fragment Analyser
   System; Agilent). The constructed library was sequenced using PE150
   sequencing mode on an Illumina Nova 6000 platform following the
   manufacturer’s instructions.

Single-nucleus data analysis

   Illumina sequencing outputs were generated for the samples as binary
   base call (BCL) files. Using the 10X Genomics Cell Ranger software
   pipeline (version 3.1.0), the BCL files were demultiplexed and
   converted into FASTQ files. The raw snRNA-seq FASTQ files were
   processed to produce a matrix of gene counts versus nuclei. Nuclei
   expressing more than 200 genes were retained. Moreover, nuclei with
   more than 5% of the expressed mitochondrial genes were filtered out.
   Library size normalization was performed on the filtered nuclei using
   the R package Seurat (version 4.0.0) to obtain the normalized count
   matrix. Principal component analysis (PCA) was performed to reduce the
   dimensionality of the log-transformed gene-barcode matrices of the top
   variable genes. The dimension reduction algorithm mutual nearest
   neighbours (MNN) was used to eliminate batch effects. Analyses of
   dimension reduction clustering and cluster-specific marker detection
   were performed with Seurat. Uniform manifold approximation and
   projection (UMAP) using the Loupe cell browser (10X Genomics) was used
   for initial visualization of the optimal cell population through
   variable gene selection. The generated clusters were identified as
   specific cell types based on the cell type-specific marker genes, which
   entailed the canonical markers for individual muscle cell types or the
   top DEGs in each cluster. The clusters were further classified into
   subclusters based on the heterogeneity of the specific cell types. DEGs
   analysis was performed for each of the cell clusters using Seurat. A P
   value < 0.05 and |log[2]foldchange|>0.58 were set as the thresholds for
   significant differential expression. Kyoto Encyclopedia of Genes and
   Genomes (KEGG) pathway enrichment analysis of the DEGs was considered
   significant at an adjusted P value < 0.05. DEGs and KEGG analyses
   revealed key candidate genes and molecular pathways in the regulatory
   network related to diaphragm fibrosis and atrophy. To infer dynamic
   gene expression and progression among the clusters and subclusters, we
   utilized the R package Monocle (version 2.9.0) to perform pseudotime
   trajectory analysis. We used the differential gene test function of the
   Monocle package to select genes (q < 0.01) that were likely to be
   informative in the ordering of cells along the pseudotime trajectory.
   We then applied protein‒protein interaction (PPI) analysis to the hub
   and differential genes, which indicated that genes situated at the
   center of the interaction network were responsible for their key
   connections with other genes.

High-dimensional WGCNA

   To investigate the intrinsic properties of the snRNA-seq data, hdWGCNA
   was performed to decipher the hub genes expression. Following the
   workflow for conducting hdWGCNA, we constructed gene coexpression
   networks based on pairwise correlations among genes using snRNA-seq
   data from four diaphragm samples. The WGCNA tree has a minimum size of
   50 genes, and similar genes are integrated into one module. Modules or
   clusters of highly correlated genes were identified, and the module
   eigengenes (MEs) describing the expression patterns of coexpression
   modules were calculated. The genes with the highest eigengene-based
   connectivity (kME) values in the modules were identified as hub genes.

Quantitative real-time polymerase chain reaction (qRT-PCR)

   To confirm the results of the snRNA-seq analysis, the expression of six
   genes, 6-Phosphofructo-2-kinase/fructose-2,6-bisphosphatase 3 (Pfkfb3),
   Platelet-derived growth factor D (Pdgfd), C-X-C Motif Chemokine
   Receptor 2 (Cxcr2), Neuronal growth regulator 1(Negr1), Semaphorin-3a
   (Sema3a) and myocyte enhancer factor-2 C (Mef2c), were determined using
   qRT-PCR. TRIzol reagent (15596026, Ambion, Carlsbad, CA, USA) was used
   for total RNA extraction from whole-diaphragm samples according to
   standard protocols. The quality of the purified RNA samples was
   evaluated with an Agilent Bioanalyzer (Agilent Technologies, Inc.,
   Santa Clara, CA, USA), and RNA samples with an RNA integrity number
   (RIN) > 8.0 were used for gene expression measurements. Reverse
   transcription was performed with the Super Script II First-Strand
   Synthesis Kit (2690A, Takara Bio USA Inc., San Jose, CA, USA), and
   qRT-PCR was performed with the SYBR Fast Universal qPCR Kit (KK4601,
   KAPA Biosystems, Woburn, MA, USA) on a CFX-Connect 96 Real-Time PCR
   detection System (Bio-Rad, Hercules, CA, USA). The primer sequences
   (5’-3’) used were as follows: Pfkfb3, forward: AAACGCCTTAAGCTCCTAGTGAT,
   reverse: GCAACGTGCCAGTCTGATGG;Pdgfd,forward: ACCTGGGTTCAAAAGAATTTA,
   reverse: GACTTGCGATCATGGTAT;Cxcr2,forward: CAACTCCCTGGTGATGCT, reverse:
   TCCTTCACAAGCGAGACC;Negr1,forward: ACCAGAGCCTTCCATTTC, reverse:
   GGTCACAGTGCCAGATTTA;Sema3a,forward: GATGAGGAACGGAGTAGGC, reverse:
   GCAAATTGGGTGGAAAGC;Mef2c,forward: CCCAACCTCTTGCCACTG, reverse:
   TTGCCATATCCGTTCCCT. The relative gene expression was quantified by the
   ΔΔCt method, with the gene Gapdh serving as the internal control.

Western blotting analysis

   A maximum of 20 mg of diaphragm tissue was homogenized on ice in 800 ul
   of lysis buffer. Total diaphragm proteins were isolated using RIPA
   buffer (R0010, Solarbio, Science & Technology Co., Ltd., Beijing,
   China) supplemented with protease inhibitor cocktails. After
   centrifugation (12,000 rpm for 15 min), the protein concentration was
   determined using a bicinchoninic acid (BCA) protein assay kit (P0010,
   Solarbio). Subsequently, 20 µg of each sample was subjected to 10%
   SDS‒PAGE and transferred to a PVDF membrane (IPVH00010, Millipore,
   Billerica, MA). After the semidry blotting procedure (50 min, 90 V),
   the membrane was incubated for 1 h at room temperature (RT) in 5% BSA
   blocking solution, followed by overnight incubation on a shaker at 4 °C
   with primary antibodies against PFKFB3 (bs-3528R, Boster Biological
   Technology Co., Ltd., Wuhan, China), PDGFD (bs-24572R, Boster), CXCR2
   (bs-1629R, Boster), NEGR1 (bs-11095R, Boster), SEMA3A (bs-10468R,
   Boster), MEF2C (bs-4130R, Boster) and β-actin (66009-1-Ig, Proteintech,
   Wuhan, China) as a loading control. The membrane was then incubated for
   1 h at room temperature on a shaker with a horseradish
   peroxidase-conjugated goat anti-rabbit secondary antibody (BL003A;
   Biosharp, Guangzhou, China). After the bands were visualized with
   tetramethylbenzidine, the blottings were scanned and quantified using
   Tanon GIS software (version 4.2; Tanon Technology, Shanghai, China).

Statistical analysis

   Statistical analysis was performed with SPSS (version 25.0; SPSS Inc.,
   Chicago, IL, USA). Differences between two groups were evaluated using
   two-tailed Student’s t-tests. All results are presented as the
   mean ± standard deviation (SD). A P value < 0.05 was considered
   statistical significance.

Results

Blood gas and arterial blood pressure

   Blood gas analysis revealed that the PaO[2] concentration remained at
   133 ± 11 mmHg for the experimental group and 125 ± 12 mmHg for the
   control group. This was also the case for the PaCO[2] concentration
   (experimental group, 35 ± 6 mmHg; control group, 33 ± 5 mmHg) and pH
   (experimental group, 7.44 ± 0.03; control group, 7.41 ± 0.02). Arterial
   blood pressure was similar in the different groups and averaged 115 ± 4
   and 127 ± 6 mmHg for the experimental group and control group,
   respectively.

Diaphragm cell clusters and subclusters

   Our snRNA-seq data revealed the intricate transcriptional dynamics of
   the diaphragm during 24 h of mechanical ventilation. Nuclei from rabbit
   diaphragms were purified and successfully constructed on the 10X
   Genomics Platform (Fig. [52]1). We profiled 14,652 nuclei
   transcriptomes from the diaphragms of 2 mechanically ventilated rabbits
   and 2 controls for analysis following normalization and quality control
   of the snRNA-seq data. Sequencing data were plotted using UMAP to
   generate functionally enriched clusters. Subsequently, the generated
   clusters were annotated according to the expression of cell
   type-specific marker genes. The fibre type-specific marker genes were
   identified as previously reported or by the top cell type-specific
   differential genes. We identified 15 clusters annotated to 10 major
   cell types in the diaphragm (Fig. [53]2a). Two top marker genes for
   each cluster were displayed on the dot plots (Fig. [54]2b). Violin
   plots showed that the marker genes Pfkfb3 and Tbc1d1 were mainly
   expressed on the clusters 1,2,8, Tek and Emcn on the clusters 9,14
   (Fig. [55]2c). Clusters 1, 2, 8 were annotated as fibroblast 1 with the
   marker genes Pfkfb3 and Tre-2/USP6, BUB2, cdc16 domain family member 1
   (Tbc1d1) (Fig. [56]2d); Cluster 5 was defined as myocyte 1 with the
   marker genes Myh7b and Mylk3 (Fig. [57]2e); Cluster 7 was defined as
   myocyte 2 with the marker genes Myh1 and Myl1 (Fig. [58]2f); Cluster 6
   was defined as fibroblast 2 containing FAP with the marker genes
   Col15a1 and Pdgfd (Fig. [59]2g); Clusters 9 and 14 were defined as
   endothelial cells marked by Tek and Emcn (Fig. [60]2h), Cluster 9 has a
   much more distinct expression profile and was defined as postcapillary
   venule (PCV) endothelial cells; Cluster 12 was defined as neurocytes
   marked by Pak5 and Ntng1 (Fig. [61]2i); Cluster 13 was divided into
   immune cells marked by Ikzf1 and Cd163 (Fig. [62]2j) and satellite
   cells marked by Pax7 and Fgfr4 (Fig. [63]2k); Cluster 15 was annotated
   as pericytes and smooth muscle cells marked by Acta2 and Rgs5
   (Fig. [64]2l). We further divided the marked cell population into cell
   subgroups. The biomarkers included global markers, which specify cell
   types, and local markers, which specify subtypes. The diaphragm has
   three different specialized fibroblast clusters: mature fibroblasts
   that are shared by all tissue fibroblasts in the steady state,
   intermediate fibroblasts between activated and quiescent states, and
   naive fibroblasts that are activated by injury. The cell population of
   fibroblasts that increased within 24 h of mechanical ventilation was
   named fibroblast 1, which consisted of two subclusters:
   myofibroblasts-myoblasts and naive fibroblasts (Fig. [65]3a, b, c).
   Endothelial cells can be separated into two subclusters: main and PCV
   endothelial cells (Fig. [66]3d, e, f). The clusters with high Pdgfd
   expression were identified as fibroblast 2, which included two
   subclusters: mature and intermediate fibroblasts (Fig. [67]3g, h, i),
   FAP were included in fibroblast 2. The small populations of immune
   cells and satellite cells in the samples were separated into immune
   cell subclusters and satellite cell subclusters (Fig. [68]3j, k, l).

Fig. 1.

   [69]Fig. 1
   [70]Open in a new tab

   Schematic diagram for single-nucleus sequencing and diaphragm fibrosis
   formation. a Schematic diagram for nuclei purification and sequencing
   from rabbit diaphragms (By Figdraw 2.0, ID: AORYU13345). b Clustering
   of snRNA-seq data through uniform manifold approximation and projection
   (UMAP) dimension reduction analysis from a combined set of four
   samples. 15 clusters annotated into 10 cell types, namely, fibroblast
   1, fibroblast 2, myocyte 1, myocyte 2, endothelial cells, neurocyte,
   immune cell and satellite cell, pericytes and smooth muscle cell. Cells
   are coloured by major cell type assignment. The proportion of
   fibroblast 1 cells increased, and the proportion of myocytes decreased
   quickly after 24 h of mechanical ventilation in the diaphragm. c
   Schematic diagram of diaphragm fibrosis during 24 h of mechanical
   ventilation (By Figdraw 2.0, ID: ITAOR9fdde). d Schematic diagram
   showing the key gene signatures and signaling pathways associated with
   the development of diaphragm fibrosis and atrophy (By Figdraw 2.0, ID:
   WTAOWa775f).

Fig. 2.

   [71]Fig. 2
   [72]Open in a new tab

   Diaphragm cell type identification. a Clustering algorithm with a
   resolution of 0.6 identified 15 clusters. b Dot plots of key marker
   gene expression levels across cell clusters. The average expression
   (dot colour intensity) and cell percentage (dot size) of the two top
   marker genes for each cluster were compared. c. Violin plots showing
   the relative expression levels of selected genes for different
   clusters. Cells are coloured according to the expression of the
   specific marker genes. UMAP of snRNA-seq data displayed the presence of
   Fibroblast 1 containing Cluster 1, 2, 8 marked by Pfkfb3 and Tbc1d1
   (d); Myocyte 1 corresponding to Cluster 5 marked by Myh7b and Mylk3
   (e); Myocyte 2 corresponding to Cluster 7 marked by Myh1 and Myl1 (f);
   Fibroblast 2 corresponding to Cluster 6 marked by Col15a1 and Pdgfrα
   (g); Endothelial cell including Cluster 9,14 marked by Tek and Emcn
   (h); Neurocyte corresponding to Cluster 12 marked by Pak5 and Ntng1
   (i); Immune cell corresponding to Cluster 13 marked by Ikzf1 and Cd163
   (j); Satellite cell corresponding to Cluster 13 marked by Pax7 and
   Fgfr4 (k); Pericytes and smooth muscle cell corresponding to Cluster 15
   marked by Acta2 and Rgs5 (l).

Fig. 3.

   [73]Fig. 3
   [74]Open in a new tab

   Subcluster and cell type analysis. Clustering and identification of
   subpopulations. Specific cell populations, such as Fibroblast 1 (a),
   Endothelial cell (d), Fibroblast 2 (g), Immune and Satellite cell (j),
   were marked with Loupe software, while other cell populations are
   marked in grey. Fibroblast 1 has two main subclusters (b):
   Myofibroblasts/myoblasts, which are marked by the genes Tns1 and Myod1,
   and naive fibroblasts, which are marked by Sec24d and Fbxo32 (c).
   Endothelial cells also had two main subclusters (e): endothelial cells
   marked by the Adamts12 and Foxp2 genes and PCV endothelial cells marked
   by Dach1 and Mecom (f). Fibroblast 2 has two main subclusters (h):
   intermediate fibroblasts, which are marked by the genes Pdk4 and
   Zbtb16, and mature fibroblasts, which are marked by Dtx3 and Phldb2
   (i). There are two main subclusters of immune and satellite cells (k):
   satellite cells, marked by the genes Fgfr4 and Pax7, and immune cells,
   marked by Ikzf1 and Cd163(l).

DEGs and KEGG analysis

   We performed DEGs and KEGG pathway enrichment analysis (Fig. [75]4). As
   shown in the volcano plots, the expression of Pfkfb3, Tbc1d1, Pdgfd,
   and Cxcr2 was significantly upregulated in the experimental group
   compared with the control group in the fibroblast 1, fibroblast 2 and
   immune cells (P < 0.05), while the expression of Ccl21, Mef2c, and
   Negr1 was significantly downregulated in the experimental group
   compared with the control group in the endothelial cells, myocyte and
   neurocyte respectively (P < 0.05). KEGG analysis revealed that the DEGs
   in the fibroblast 1, fibroblast 2, immune cell, endothelial cell,
   myocyte and neurocyte were functionally enriched in the insulin
   signaling pathway, phosphatidylinositol 3-kinase (PI3K) -Akt signaling
   pathway, phosphalipase D (PLD) and Rap1 signaling pathway, chemical
   carcinogenesis-reactive oxygen species, calcium signaling pathway and
   leukocyte transendothelial migration.

Fig. 4.

   [76]Fig. 4
   [77]Open in a new tab

   Differentially expressed genes and Kyoto Encyclopedia of Genes and
   Genomes analysis. Volcano plots showing the top twenty up- or
   downregulated genes for each cell type. The red dots represent
   upregulated genes, and the blue dots represent downregulated genes. (P
   value < 0.05, and |log2foldchange| > 0.58). Bubble plot indicating the
   top enriched pathways for each cell type based on KEGGpathway
   enrichment analysis of differentially expressed
   genes([78]www.kegg.jp/kegg/kegg1.html). The sizes of the dots represent
   the number of genes included in each pathway. The colour gradient of
   dots represents the adjusted P values of each enriched pathway. The
   genes Pfkfb3, Tbc1d1 and the insulin signaling pathway; Pdgfd and the
   PI3K-Akt signaling pathway; Cxcr2 and PLD, Rap1 signaling pathways;
   Ccl21 and chemical carcinogenesis-reactive oxygen species signaling
   pathway; Mef2c and the calcium signaling pathway; Negr1 and leukocyte
   transendothelial migration signaling pathway are labelled by the black
   and red bars in Panels a, b, c, d, e, and f.

hdWGCNA and PPI network

   We performed hdWGCNA on a single-nucleus transcriptomic dataset of
   diaphragm tissue from two mechanically ventilated rabbits and two
   control ones, which revealed distinct network structures and sets of
   gene modules. After weighting the median connectivity and scale-free
   topology model fit, a power value of 6 was selected. The final analysis
   resulted in eight modules that were derived based on a scale free
   net-work structure. Within cell type-specific modules, the hub genes
   exhibited the highest intramodular connectivity among the coexpression
   networks and displayed biologically relevant information for each major
   cell type. We identified brown module was consistent with the
   fibroblast 1 cluster, containing two labelled hub genes, Pfkfb3 and
   Tbc1d1; the green module was consistent with the fibroblast 2 cluster,
   containing a labelled hub gene, Pdgfd; the turquoise module was
   consistent with the myofiber cluster, containing a labelled hub gene,
   Ryr3; and the yellow module was consistent with the endothelial cell
   cluster, containing a labelled hub gene, Sema3a (Fig. [79]5). Our
   analysis ultimately led to the identification of five hub genes:
   Pfkfb3, Tbc1d1, Pdgfd, Ryr3, and Sema3a. We constructed nine PPI
   networks centred on the proteins PFKFB3, TBC1D1, PDGFD, CXCR2, SEMA3A,
   CCL21, NEGR1, MEF2C, and RYR3 to determine their relationships with
   other proteins (Fig. [80]6).

Fig. 5.

   [81]Fig. 5
   [82]Open in a new tab

   High-dimensional weighted gene coexpression network analysis (hdWGCNA).
   a The top left panel depicts the soft power threshold for choosing a
   scale-free topology model. The average connectivity of the topological
   network was most stable at the lowest soft threshold of 6. b Dynamic
   Tree Cut algorithm for gene clustering. Each leaf on the tree
   represents a gene, and the colour at the bottom indicates the
   assignment to a specific coexpression module. c Bubble plots displaying
   the scores obtained for eight modules in eight cell subtypes. d
   snRNA-seq UMAP coloured by module eigengene (ME) for eight coexpression
   modules. e The eight modules’ values of -log10 (Adj P value) on the
   Y-axis versus average log2 (fold change) on the X-axis are shown for
   each cell type, including endothelial cells, fibroblast 1, fibroblast
   2, myocytes, neurocytes, pericytes, satellite cells and smooth muscle
   cells. The larger the value is, the greater the fit between the
   coloured module and the specific cell type. The plot shows that
   fibroblast 1 features a brown module, fibroblast 2 features a green
   module, endothelial cells feature a yellow module, and myocyte cells
   feature a turquoise module. f UMAP plot of the gene coexpression
   network with eight coloured modules and the top two hub genes labelled
   in each module. Nodes are coloured according to coexpression module
   assignment. g Hub genes in each module were identified by
   eigengene-based connectivity (KME). The following five hub genes
   associated with fibrosis were identified: Pfkfb3 and Tbc1d1 in the
   brown module, Pdgfd in the green module, Sema3a in the yellow module
   and Ryr3 in the turquoise module.

Fig. 6.

   [83]Fig. 6
   [84]Open in a new tab

   Protein-protein interaction network (PPI network).We constructed nine
   PPI networks displaying protein-protein interactions among the related
   genes. The nodes represent proteins, and the edges represent the
   interaction strength between two proteins. The proteins PFKFB3 (a),
   PDGFD (b), CXCR2 (c), CCL21 (d), SEMA3A (e), RYR3 (f), MEF2C (g),
   NEGR1(h) and TBC1D1 (i) which are located at the hub of the interaction
   network, are responsible for diaphragm fibrosis and atrophy.

Pseudotime trajectory analysis

   We identified four trajectories corresponding to distinct cell fates,
   termed the satellite cell, FAP cell, endothelial cell and immune cell
   trajectories. Pseudotime analysis revealed that the satellite cell
   trajectory started from Cluster 13 to Cluster 1, indicating that the
   quiescent satellite cells expressing Pax7 differentiated from
   progenitors to myoblasts expressed Myod1 after activation during
   mechanical ventilation. Monocle analysis clarified the satellite cell
   trajectory by reconstructing two main branches. One trajectory
   identified cells progressing towards myofiber differentiation, as
   indicated by the activation of Myod1. Another branch continued to
   proliferate, retaining the identity of the satellite cell (Fig. [85]7a,
   b, c). The ordering of cells in the trajectory tree suggested that the
   FAP cell trajectory seemed to start principally from Cluster 12, then
   move towards Cluster 11, and finally to Clusters 1, 2, 8 representing
   the trajectory from FAP cell (contained in fibroblast 2) to naive
   fibroblasts (contained in fibroblast 1) (Fig. [86]7d, e, f). The
   trajectory of endothelial cells started principally from Clusters 9 and
   14, then moved towards Cluster 10, and finally to Clusters 1, 2, 8
   representing the trajectory from endothelial cells to naive fibroblasts
   (Fig. [87]7g, h, i). The immune cell trajectory started principally
   from Cluster 13 and then moved to Clusters 1, 2, 8 representing the
   trajectory from immune cells to naive fibroblasts (Fig. [88]7j, k, l).

Fig. 7.

   [89]Fig. 7
   [90]Open in a new tab

   Single-nucleus trajectory analysis. Pseudotime trajectories
   reconstructed by Monocle 2 highlight the utility of nucleus-level
   resolution, which reflects the transcriptional dynamics that regulate
   the formation of diaphragm fibrosis. a UMAP-coloured only the satellite
   cell and myoblast; b the trajectory of satellite cells to myoblasts; d
   UMAP-coloured fibroblast 2, other cells and fibroblast 1; e The
   trajectory of FAP (contained in fibroblast 2) to naive fibroblasts
   (contained in fibroblast 1); g UMAP-coloured endothelial cells, other
   cells and fibroblast 1; h The trajectory of endothelial cells to naive
   fibroblasts (contained in fibroblast 1); j UMAP-coloured immune cells
   and fibroblast 1; k The trajectory of immune cells to naive fibroblasts
   (contained in fibroblast 1). The beginning and end of the timeline in
   pseudotime trajectories are displayed with the direction of the arrows
   (c, f, i, l).

Quantitative real-time PCR and Western blotting analysis

   The qRT-PCR data showed that the mRNA levels of the Pfkfb3, Pdgfd,
   Cxcr2, and Negr1 genes were significantly higher in the mechanically
   ventilated group than in the control group (P < 0.01), and the
   expression levels of Sema3a and Mef2c were significantly lower
   (P < 0.01). Western blotting analysis revealed that the protein
   expression of PFKFB3, PDGFD, and CXCR2 tended to increase and that the
   protein expression of SEMA3A, NEGR1 and MEF2C tended to decrease in the
   experimental group compared with the control group, but these
   differences did not reach statistical significance (Fig. [91]8).

Fig. 8.

   [92]Fig. 8
   [93]Open in a new tab

   Quantitative real-time polymerase chain reaction (qRT-PCR) and Western
   blotting. qRT-PCR and Western blotting showing the expression of
   selected genes and corresponding proteins in the mechanically
   ventilated diaphragm and control ones. The bar graphs show the
   quantification of the mRNA expression of Pfkfb3 (a), Pdgfd (b), Cxcr2
   (c), Negr1 (d), Sema3a (e), and Mef2c (f) normalized to that of GAPDH.
   Asterisks indicate significant differences (n = 3 in each group)
   (*P < 0.05). Western blotting analysis of the protein expression of
   PFKFB3 (g), PDGFD (h), CXCR2 (i), NEGR1 (j), SEMA3A (k), and MEF2C (l)
   normalized to that of β-actin as a loading control.

Discussion

   This study is the first attempt to gain insights into the cellular and
   molecular mechanisms underlying the early changes that occur in the
   diaphragm in response to mechanical ventilation by snRNA-seq using a
   rabbit model. We found that after 24 h of mechanical ventilation
   diaphragm cell changes included a decreased proportion of myofibres and
   an increased proportion of fibroblasts, which may contributed to the
   development of diaphragm fibrosis. In contrast to the traditional
   perspective, the study found evidence of myocyte atrophy and fibrosis
   in the diaphragm, as well as potential reductions in diaphragmatic
   contractile force. These findings suggest that myocyte atrophy and
   fibrosis may jointly contribute to VIDD, but additional studies are
   needed to confirm their causality and to investigate other potential
   factors involved in this condition.

   During 24 h of mechanical ventilation in rabbits, the activation of
   fibroblasts and their subsequent differentiation into myofibroblasts,
   which depends on a combination of specific factors induced the
   development of fibrosis. The genes Pfkfb3 and Tbc1d1 and the insulin
   signaling pathway may drive glycolysis to meet the energy requirements
   for fibroblast proliferation during mechanical ventilation. In
   addition, proliferation of FAP and fibroblast-to-myofibroblast
   transdifferentiation (FMT) associated with the corresponding genes
   Pdgfd and the PI3K-Akt signaling pathway, EndMT associated with the
   corresponding gene Sema3a and the signaling pathway chemical
   carcinogenesis-reactive oxygen species, immune cell infiltration
   associated with the regulatory gene Cxcr2, PLD and Rap1 signaling
   pathways are the three main processes involved in early stage of
   driving diaphragm fibrosis. Phrenic nerve ending loss resulting in a
   lack of trophic support with decreased Negr1 gene expression and
   leukocyte transendothelial migration signaling pathway activation,
   vascular reduction leading to diaphragm hypoperfusion with decreased
   Ccl21 gene expression, calcium signaling pathway activation with
   decreased Mef2c and Ryr3 genes expression are the three main factors
   responsible for myocyte atrophy. Muscle satellite cells are progenitors
   of muscle fibres and reside under the basal lamina in direct contact
   with myofibers^[94]19. In rabbits, the diaphragm muscle shows earlier
   and more severe fibre damage during mechanical ventilation. Upon
   injury, quiescent satellite cells began to proliferate and
   differentiate into myofibroblasts-myoblasts with upregulated Myod1
   expression, which was confirmed by the satellite cell trajectory
   analysis in our study. However, relatively few myogenic precursors
   expressing Myod1 exhibited a limited proliferation and differentiation
   capacity to meet functional demands as diaphragmatic fibres are lost.
   Under these conditions, normal repair of the diaphragm is replaced by
   fibrosis, which is a pathological repair process.

   Aerobic glycolysis has been increasingly recognized as an important
   pathogenic process that contributes to the development of many types of
   fibrosis^[95]20. Metabolic reprogramming of glycolysis is necessary to
   meet the energetic demand of activated fibroblasts compared to
   quiescent fibroblasts, and this metabolic switch highlights the
   importance of glycolysis for fibrosis development^[96]21. Pfkfb3 is a
   critical rate-limiting enzyme that promotes glycolysis metabolism by
   modulating and maintaining the intracellular concentration of
   fructose-2,6-bisphosphate (F-2,6-BP)^[97]22. Upregulation of Pfkfb3
   increased glycolysis through the TGF-β signaling pathway, and the
   production of lactate directly contributed to fibroblast activation and
   proliferation in renal fibrosis models of unilateral ureter obstruction
   and kidney ischaemia/reperfusion^[98]23. The overexpression of Pfkfb3,
   which promotes glycolysis in cardiac fibroblasts, results in cardiac
   fibrosis while postmyocardial infarction^[99]24. In an experimental
   model of pulmonary fibrosis, glycolytic reprogramming promoted lung
   myofibroblast activation by upregulating the expression of the critical
   glycolytic enzyme Pfkfb3^[100]25. Tbc1d1 is a member of the TBC1
   Rab-GAP family of proteins and has distinctive roles in regulating both
   insulin- and contraction-stimulated glucose transport in skeletal
   muscle^[101]26. Insulin and exercise stimulate the phosphorylation of
   Tbc1d1 by upregulating insulin receptor substrate (IRS)-PI3K and
   increasing AMP-activated protein kinase activity, and Tbc1d1
   phosphorylation facilitates the translocation of glucose transporter 4
   (GLUT4) to the cell surface membrane, resulting in increased cellular
   glucose uptake^[102]27. Pfkfb3 and Tbc1d1 are the critical nodes
   involved in the insulin signaling pathway and are central regulators of
   metabolic function in skeletal muscle. The increased expression of
   Pfkfb3 and Tbc1d1 in our snRNA-seq data suggested that there were
   strong connections between the transcriptomic and metabolomic profiles
   of the diaphragm in mechanically ventilated rabbits for 24 h. We
   speculated that enhanced glucose metabolism driven by Pfkfb3, Tbc1d1
   and the insulin signaling pathway may produce a large amount of energy
   to meet the demand of rapidly differentiating fibroblasts, which is
   necessary for the development of diaphragm fibrosis.

   FAP are nonmyogenic mesenchymal stem cells that reside in the
   interstitial space between myofibres and are characterized by the
   expression of Pdgfrα, which induces FMT and increases the production of
   ECM proteins for muscle structure support^[103]28,[104]29. Pdgfd
   strongly activated the PI3K-AKT pathway. AKT-catalysed phosphorylation
   of the forkhead box gene, group O1 (Foxo1), results in its
   translocation from the nucleus to the cytosol and prevents
   Foxo1-dependent gene expression of manganese superoxide dismutase
   (Mnsod). The downregulation of Mnsod leads to the excess production of
   ROS and FMT^[105]30. Previous studies reported that FAP, which are the
   predominant contributors to pathological diaphragm fibrosis, expressed
   greater levels of the Pdgfrα protein in the diaphragm than in the limb
   muscles of mdx mice (a muscular dystrophy model)^[106]31. In the
   present study, increased Pdgfd expression, as determined by DEGs
   analysis of single-nucleus sequencing data, was confirmed by qRT-PCR
   and Western blotting. We postulated that increasing Pdgfd expression
   together with PI3K/AKT/FoxO1 pathway activation were the main drivers
   of diaphragm fibrosis in response to muscle injury during mechanical
   ventilation based on the DEGs analysis via snRNA-seq.

   Mechanical ventilation triggered the recruitment of circulating immune
   cells into diaphragmatic tissue after diaphragm injury. In our study,
   the increase in the expression of Cxcr2 and activation of the Rap1 and
   PLD signaling pathways could clarify the molecular mechanisms by which
   immune cells infiltrate diaphragmatic tissue and transform into
   fibroblasts. Cxcr2 belongs to the γ subfamily of GTP binding proteins
   (G-proteins) and is the receptor for CXCL8. Once the chemokine CXCL8
   binds to the receptor Cxcr2, the β and γ subunits of the G-protein are
   released and the downstream small GTPase RAS-related protein 1 (Rap1)
   is activated, PLD is also activated to produce the second messenger
   phosphatidic acid (PA), which increases immune cell adhesion and
   migration, this two processes are closely related to the recruitment of
   immune cells in response to acute diaphragm injury^[107]32,[108]33.

   Vascular endothelial cells form a monolayer lining the vasculature
   throughout the body and play diverse roles in maintaining internal
   homeostasis. EndMT is a transcriptional program that downregulates the
   expression of endothelial cell marker genes and upregulates the
   expression of mesenchymal cell marker genes. Endothelial cells can
   transdifferentiate into fibroblasts via EndMT and contribute to
   extracellular matrix deposition by driving fibrogenesis^[109]34. SEMA3A
   is a special membrane-associated secreted protein with various
   biological properties, and the SEMA3A/LIMK/p-cofilin/actin axis is
   involved in EndMT. Endothelial cells acquire migratory capabilities
   after actin remodelling during EndMT. SEMA3A blocks TGF-β-induced EndMT
   by inhibiting cytoskeletal actin remodelling and cell migration
   mediated by LIMK/p-cofilin signaling^[110]35. In our study, reduced
   expression of Sema3a/SEMA3A, as determined by qRT-PCR and Western
   blotting, increased LIMK/p-cofilin signaling and contributed to
   TGF-β-induced EndMT via dynamic actin remodelling, thus contributing to
   the development of diaphragm fibrosis.

   Taken together, these findings indicate that FAP are the most important
   contributors to mechanical ventilator-triggered fibroblast
   transformation and activation as early as 24 h. The contribution of
   immune cells and endothelial cells to diaphragm fibrosis is likely
   limited due to the rarity of immune cells, and EndMT is a complex
   phenotypic transformation process.

   Negr1 is a glycosylphosphatidylinositol (GPI)-anchored membrane
   synaptic adhesion molecule highly expressed in the cerebral cortex and
   hippocampus^[111]36. Negr1 promotes sympathetic axonal growth and
   affects adult neurogenesis by regulating leukocyte inhibitory factor
   receptor (LIFR) signaling^[112]37. In addition, leukocyte
   transendothelial migration (LTM) is an inflammatory process in which
   leukocytes cross through the blood vessel wall, leading to tissue
   damage and organ failure^[113]38. LTM induces the apoptosis of neurons
   by secreting various cytotoxic molecules, such as C-X-C motif chemokine
   ligand 10 (Cxcl10)^[114]39. In our study, reduced secretion of Negr1
   combined with LTM activation resulted in the progressive loss of
   phrenic nerve terminals through inflammatory cell recruitment and
   apoptosis promotion during 24 h of mechanical ventilation in the
   diaphragm. The degeneration of phrenic neuron endings contributes to a
   lack of trophic support for muscle tissue and diaphragm dysfunction.

   Angiogenesis refers to the formation of new blood vessels that depend
   on endothelial cell activation and proliferation via the release of
   endothelial growth factors^[115]40. Angiogenesis plays a key role in
   wound healing during the inflammation phase by facilitating the
   recruitment of inflammatory cells and the secretion of growth
   factors^[116]41. Ccl21 is the main chemokine that promotes angiogenesis
   by enhancing the formation of blood vessels through binding to CCR7 and
   activating the PI3K pathway^[117]42. In our study, Ccl21 expression was
   reduced at the 24-hour time point in the experimental group compared to
   the control group. Correspondingly, we propose that reduced Ccl21
   expression results in blood vessel loss, especially blood capillary
   loss, and decreased blood flow leading to diaphragm dysfunction within
   24 h of mechanical ventilation.

   Previous animal studies revealed that slow and fast-twitch diaphragm
   fibres exhibit atrophy simultaneously induced by MV^[118]43, while our
   study revealed that slow fibre atrophy is induced by MV via the calcium
   signaling pathway and the genes Ryr3 and Mef2c. Defective intracellular
   Ca^2+ handling is central to depressed diaphragm contractility in
   mechanically ventilated rabbits^[119]44,[120]45. There are three
   different isoforms of the ryanodine receptor Ca^2+ release channel
   (Ryr1, Ryr2, and Ryr3) expressed on the sarcoplasmic reticulum in
   skeletal muscle, and the Ryr1 isoform represents the main isoform of
   the Ca^2+ release channel^[121]46. In VIDD, diaphragmatic mitochondrial
   oxidative stress causes Ryr1 channel complex instability and Ca^2+
   leakage, then cytoplasmic Ca^2+ overload leads to calpain activation
   and endoplasmic reticulum (ER) stress, which contribute to myocyte
   proteolysis and apoptosis^[122]47,[123]48. Ryr3 expression is found
   only in the slow fibres of the diaphragm in adult rabbits^[124]49. One
   study showed that RyR1 and RyR3 coexpressed in the same cell in rabbit
   diaphragm forms a homotetramer^[125]50. The coexpression of the RyR1
   and RyR3 calcium release channels in skeletal muscles may result in a
   diversification of the excitation-contraction coupling apparatus to
   obtain a fine regulation of muscle contraction^[126]51. In our study,
   the identification of the turquoise module, which includes the hub gene
   Ryr3 is a novel finding, given the relatively low expression levels of
   RyR3 compared to RyR1 in the diaphragm. We propose RyR3 might play the
   important role leading to diaphragm slow-twitch fibre atrophy and
   contractile dysfunction in the diaphragm. The myocyte enhancer factor-2
   (Mef2) family is a family of transcription factors that regulate
   skeletal muscle repair and regeneration through the activation of
   satellite cells. Motor nerve activity is mediated by
   calcineurin-dependent genes via Mef2 transcription to control
   fundamental cellular processes, thereby establishing special
   physiological myofiber subtypes^[127]52. In vertebrates, this family
   has four members: Mef2a, Mef2b, Mef2c and Mef2d. Mef2c, an essential
   transcription factor, regulates the transcription of slow-twitch fibre
   specific genes, particularly Myh7, and its product MyHC-β/slow^[128]53.
   Prolonged bed rest in humans induces robust inhibition of the Mef2
   family members, particularly Mef2c, which ultimately reduces skeletal
   muscle mass^[129]54. In our study, phrenic nerve ending loss led to a
   decrease in nerve impulses combined with mechanical ventilator
   unloading in the diaphragm, which induced Mef2c downregulation and
   slow-twitch fibre atrophy.

   In contrast to the traditional perspective, we speculated that
   diaphragm atrophy may originated from three main important factors,
   namely, phrenic nerve ending loss and trophic support deletion with the
   corresponding gene Negr1, diaphragm blood vessel reduction and
   hypoperfusion with the corresponding gene Ccl21, and calcium signaling
   pathway activation with decreased Mef2c and Ryr3 expression.

Limitations

   There are several limitations in our study. First, the data emphasized
   that mechanical ventilation contributed to acute diaphragm fibrosis as
   early as after 24 h. However, fibrosis is a chronic lethal disorder and
   requires long-term progression. The extension of research time may help
   identify additional information that underlies fibrosis progression.
   Second, snRNA-seq in multinucleated cells is complex because we cannot
   determine which nuclei originate from the same cell. Additionally,
   isolating the nucleus means that only a fraction of the transcripts
   within a cell are surveyed; nuclear resident transcripts such as long
   noncoding RNAs (lncRNAs) are not included in this type of analysis.
   Third, there were opposite trends in Negr1 mRNA and NEGR1 protein
   expression. We speculated that the increase in Negr1 mRNA expression
   may reflect a transient response followed by decreased expression due
   to posttranslational modifications. Finally, we did not assess
   diaphragm muscle fiber injury, atrophy, or diaphragm strength.
   Therefore, our findings can only be correlated with 24-hour mechanical
   ventilation in rabbits and not directly with VIDD, atrophy, or injury.

Conclusions

   Our study found that 24 h of mechanical ventilation in rabbits led to
   changes in FAP proliferation, EndMT, and immune cell infiltration in
   the diaphragm. However, it is important to note that our findings are
   speculative and cannot be definitively linked to VIDD due to the
   limitations of our study. Specifically, we did not perform functional
   evaluations to confirm the presence of VIDD, nor did we have data to
   validate diaphragm unloading during the ventilation period. Despite
   these limitations, our findings suggest that FAP proliferation and
   immune cell infiltration may play a role in the early stage of driving
   diaphragm fibrosis during mechanical ventilation. Future studies are
   needed to investigate the potential mechanisms underlying these changes
   and their potential contribution to VIDD.

Electronic supplementary material

   Below is the link to the electronic supplementary material.
   [130]Supplementary Information 1.^ (81.8KB, jpg)
   [131]Supplementary Information 2.^ (56.6KB, jpg)
   [132]Supplementary Information 3.^ (107.1KB, jpg)
   [133]Supplementary Information 4.^ (72.8KB, jpg)
   [134]Supplementary Information 5.^ (78.7KB, jpg)
   [135]Supplementary Information 6.^ (77.6KB, jpg)

Acknowledgements