RNA-sequencing demonstrates transcriptional differences between human vocal fold fibroblasts and myofibroblasts

Abstract

Background

Differentiation of fibroblasts into myofibroblasts is necessary for wound healing, but excessive myofibroblast presence and persistence can result in scarring. Treatment for scarring is limited largely due to a lack of comprehensive understanding of how fibroblasts and myofibroblasts differ at the transcript level. The purpose of this study was to characterize transcriptional profiles of injured fibroblasts relative to normal fibroblasts, utilizing fibroblasts from the vocal fold as a model.

Results

Utilizing bulk RNA sequencing technology, we identified differentially expressed genes between four cell lines of normal fibroblasts (cVFF), one line of scarred fibroblasts (sVFF), and four lines of fibroblasts treated with transforming growth factor-beta 1 (TGF-β1), representing an induced-scar phenotype (tVFF). Principal component analysis revealed clustering of normal fibroblasts separate from the clustering of fibroblasts treated with TGF-β1; scarred fibroblasts were more similar to normal fibroblasts than fibroblasts treated with TGF-β1. Enrichment analyses revealed pathways related to cell signaling, receptor-ligand activity, and regulation of cell functions in scarred fibroblasts, pathways related to cell adhesion in normal fibroblasts, and pathways related to ECM binding in fibroblasts treated with TGF-β1. Although transcriptomic profiles between scarred fibroblasts and fibroblasts treated with TGF-β1 were relatively dissimilar, the most highly co-expressed genes were enriched in pathways related to actin cytoskeleton binding, which supports the use of fibroblasts treated with TGF-β1 to represent a scarred cell phenotype.

Conclusions

Transcriptomics of normal fibroblasts differ from myofibroblasts, including from those retrieved from scar and those treated with TGF-β1. Despite large differences in transcriptomics between tVFF and sVFF, tVFF serve as a useful in vitro model of myofibroblasts and highlight key similarities to myofibroblasts extracted from scar pathology, as well as expected differences related to normal fibroblasts from healthy vocal folds.

Background

Wound healing is a complex process that occurs throughout the human body. When an injury occurs, fibroblasts migrate to the injury site and become activated, resulting in differentiation into specialized cells called myofibroblasts. Myofibroblast differentiation is a physiological response required for wound healing, as these cells produce increased extracellular matrix (ECM) proteins – particularly fibrous collagen – to achieve contracture and wound closure [1]. Myofibroblasts are identified by increased expression of α-smooth muscle actin (αSMA), providing specialized stress fibers that allow for the contractility required in wound healing [2]. Without myofibroblast presence and response, poor healing and chronic wounding becomes a concern. However, excess persistence and presence of these activated cells can also result in a pathological response, known as scar or fibrosis. The role and balance of myofibroblast actions can be delicate – biological signals are meant to trigger apoptosis following closure of the wound, but any persistence of myofibroblast activity will lead to excess contracture [3]. Scarring is challenging to treat and leaves devastating and pervasive effects on function, which have been observed throughout the body, including in the skin, liver, lungs, heart, and larynx, among others.

Fibroblasts found in the vocal folds experience a similar process of wound healing following injury, and there is a lack of consistent and successful treatment for vocal fold scarring. Fibroblasts are found in the vocal fold lamina propria, which is a layer of tissue just deep to the epithelium. Unrestricted mucosal wave of the epithelium and lamina propria (termed the “cover” layer of the vocal folds) is what allows a voice to be clear and dynamic, without hoarseness or perceived difficulty in voice production. If vocal fold scarring occurs in this cover layer, vocal fold vibration is disrupted, and effects on voice quality can be detrimental. Vocal fold fibroblasts (VFF) are exposed to significant mechanical stress through vibration that occurs with all human vocalizations. Given that fibroblasts are crucial in regulating ECM production and degradation, these cells are responsible for maintaining the mechanical properties and structure throughout this mechanical stress exposure [4]. With excessive mechanical stress, as well as the position of the vocal folds within the airway causing exposure to environmental irritants, changes to ECM regulation can occur, which lead to pathogenesis in the vocal fold lamina propria, such as vocal fold scarring [5]. Because of this unique positioning in the airway and exposure to constant, and possibly harmful, mechanical loading, VFF can serve as a model for fibroblast response to injury and scarring that also occurs in other areas of the body.

One challenge related to establishing successful treatments for scarring results from a lack of comprehensive understanding of how quiescent fibroblasts and activated myofibroblasts differ at the transcript level. For example, the signals that trigger apoptosis of myofibroblasts following wound closure and the timing of these signals are unknown [6]. In order to investigate differences in transcription between cells from uninjured tissues and those from scarred tissues, high-throughput methods to assess the entire transcriptome, such as RNA sequencing (RNAseq), have been employed in the last decade as a result of improving technology. RNAseq has most notably been utilized in skin research [7,8,9], likely due to the ease of accessibility to the skin and aberrant healing patterns that may not be accessible in other areas of the body. Yet, there is a paucity of research comparing the transcriptomic profile of normal, uninjured fibroblasts to activated fibroblasts (i.e., myofibroblasts) found in scarring, in other areas of the body, particularly the vocal folds. One recent study investigated bulk RNAseq of human VFF treated with transforming growth factor beta 1 (TGF-β1) in order to identify pathways that may play a role in myofibroblast activation; this research focused primarily on the role of the myocardin related transcription factor a (MRTF-A)/serum response factor (SRF) axis in myofibroblast activation and function [10].

The purpose of this study was to comprehensively characterize gene expression of injured fibroblasts (differentiated into myofibroblasts) relative to normal fibroblasts. Given the challenge of acquiring pathologic fibroblasts from vocal fold scarring, characterization of gene expression of an in vitro model of scarred fibroblasts is also required to allow for therapeutic investigation in the absence of pathological cells. It has been well-established in prior literature that treating normal vocal fold fibroblasts with TGF-β1 in cell culture conditions leads to differentiation of those fibroblasts into myofibroblasts [11, 12]. This “forced” differentiation provides an avenue to study myofibroblast gene expression and function without requiring cells retrieved from pathology.

In the present study, transcriptomics of one rare cell line of human patient-derived scarred vocal fold fibroblasts (sVFF) relative to four lines of control (healthy) VFF (cVFF) were analyzed utilizing bulk RNAseq. Additionally, these groups were compared to the same four lines of human VFF following TGF-β1 treatment (tVFF), representing a scar model. These comparisons allow the identification of key differentially expressed (DE) genes between groups and will elucidate potential gene transcripts involved in myofibroblast differentiation and scar pathophysiology. Results will provide a comprehensive understanding of DE genes and potential biological pathways involved in scarring.

Methods

Cell culture

Five primary human VFF lines were utilized – four cVFF lines from patients without any evidence of laryngeal pathology or disease and one sVFF line from a patient who was already undergoing VF surgery for vocal fold scarring. The cVFF lines were from a 21-year-old male, 41-year-old female, 58-year-old male, and 59-year-old female. They were obtained from procurement of autopsy samples from healthy cadavers at University of Wisconsin-Madison, Department of Pathology and Laboratory Medicine, as previously reported and characterized [4]. The sVFF line was from a 56-year-old female who was previously diagnosed with VF scarring and undergoing laryngofissure with bilateral mucosal grafts as management for her severely scarred VFs. This sVFF line was also previously characterized [13].

VFF were recovered from cryopreservation at either passage 2 or 3, with passaging ratios of 1:3 utilized for seeding of cells into 10 cm culture dishes at defined cell densities based upon manual counting of trypsinized cell suspensions (0.25% Trypsin-EDTA, Gibco, Grand Island, NY) using a hemocytometer and subsequent cell concentration calculations. These cells were subcultured in media containing Dulbecco’s modified Eagle’s medium (DMEM; Cytiva, Logan, UT) supplemented with 10% fetal bovine serum (FBS; Peak Serum, Fort Collins, CO), 1% penicillin-streptomycin (Cytiva), and 1% non-essential amino acids (Sigma). After the desired concentration of cells was achieved, cells were plated at the appropriate density. All cells were at passage 5 when plated for the experiment. Each of the four cVFF lines was plated across multiple culture dishes to allow for our experimental groups: cVFF (untreated), tVFF (scar model), and nVFF (neutralization). sVFF was also plated and split into an untreated sVFF group and neutralized sVFF group. Neutralized sVFF were subjected to neutralization with the TGF-β-1,2,3 blocking antibody (Invitrogen, Carlsbad, CA) to investigate whether there were any significant effects from TGF-β neutralization in already scarred cells while neutralized healthy VFF could investigate whether transcriptomic changes in tVFF were specifically due to the TGF-β1 treatment.

One day after plating cells, cells were starved for 24 h in FBS-free media to arrest additional growth. After serum starvation, cells received the appropriate media treatment based on experimental group. cVFF and sVFF received media containing 2.5% FBS (Sigma) without additional treatment. tVFF received media containing 2.5% FBS and 10 ng/ml of TGF-β1 (Sigma, St. Louis, MO); this concentration was selected based on prior literature demonstrating increased αSMA expression [11], increased collagen secretion, and increased fibroblast migration in culture [14]. nVFF and neutralized sVFF received media containing 2.5% FBS and 1 µg/ml of TGFβ 1,2,3 blocking antibody (Invitrogen). All cells were maintained in treatment media conditions for 5 days (day 2 to day 7), with media changing on day 4, as previously described [12] and displayed in Fig. 1.

Fig. 1

Experimental groups and timeline of fibroblast cultures. At day 0, previously expanded healthy VFF and scar VFF were split and plated in experimental groups (passage 5). Following plating, day 1 required 24 h of serum starvation prior to any treatment. On day 2, new media containing low serum and applicable treatment was added to the cell. Untreated cells, including cVFF (green) and sVFF (purple) received low serum media only; TGF-β1 treated cells, tVFF (blue), received low serum media + 10 ng/ml TGF-β1; and TGF-β neutralized cells, nVFF (pink), received low serum media + 1 µg/ml TGF-β-1,2,3 blocking antibody. On day 4, cells underwent media change with continued applicable treatment. On day 7, cells were harvested for RNA extraction. Created with BioRender.com

Immunofluorescence staining

Immunohistochemical staining for αSMA was performed to confirm that tVFF had differentiated into myofibroblasts while cVFF were not differentiated. Following the treatment with or without the addition of TGF-β1 or neutralizing antibody, coverslips from 12-well plates were washed with PBS and fixed in 4% paraformaldehyde for 10 min at room temperature (RT). Coverslips were subsequently washed with a wash buffer containing 0.1% Triton X-100 (Sigma) and then underwent permeabilization with a 0.5% Triton X-100 buffer for 10 min at RT. Following permeabilization, a blocking solution made from 10% goat serum and 2% bovine serum albumin (BSA; ThermoFisher) in PBS was applied for 30 min at RT. Coverslips were then incubated with mouse αSMA primary antibody (Sigma) at 1:250 dilution with 2% BSA in PBS, for 90 min, at RT. After 90 min, coverslips were washed three times with the wash buffer of 0.1% Triton X-100. Next, secondary anti-mouse antibody (goat A488, Invitrogen) diluted to 1:200 in 2% BSA was applied to the coverslips and incubated for 60 min at RT and in the dark, to preserve fluorescence. Following, the 60-min incubation, coverslips were again washed three times with wash buffer. Slides were then prepared with one drop of Vectashield with DAPI (Vector Laboratories, Newark, CA) and coverslips were flipped over and mounted on slides. Final fluorescent images were captured on a Nikon Eclipse E600 fluorescence inverted microscope with an Olympus DP 71 camera and associated DP71 software, Olympus Corporation.

Total RNA isolation

To obtain total RNA, VFF lines were trypsinized and harvested at a density no less than 5.0 × 10⁵ cells/dish, to ensure adequate final RNA concentration for downstream RNAseq. Cells were then centrifuged at 1000 rpm for 10 min to retrieve the cell pellet. Total RNA was purified from cell pellet samples using the RNeasy Plus Mini Kit (Qiagen, Valencia, CA) according to the manufacturer’s protocol. Initial quantification of total RNA samples was performed using a NanoDrop 1000 spectrophotometer (ThermoFisher Scientific, Waltham, MA), and integrity of samples was confirmed with the following three criteria for inclusion for RNAseq: concentration > 50 ng/ml, A260 and A280 ratings between 1.8 and 2.1, and A260/A230 ratio > 1.8.

Construction and sequencing of libraries, RNA sequencing, and read mapping

Total RNA was submitted to the University of Wisconsin-Madison Biotechnology Center for verification of purity and integrity via the NanoDrop One Spectrophotometer and Agilent 4200Tapestation, respectively. Following this verification, samples were prepared using the TruSeq^® Stranded Library Preparation (Illumina Inc., San Diego, California, USA). For each library preparation, 300 ng of total RNA input was used and underwent Poly-A enrichment. The RNA was fragmented for 5 min, and random primers were added for first strand cDNA synthesis followed by second strand synthesis. Double-stranded cDNA was purified using 1.8X Agencourt AMPure XP beads. Products underwent adenylation of 3’ ends followed by ligation to Illumina undiluted unique dual indexes to increase sample throughput. These products were again purified using 0.8X Agencourt AMPure XP beads twice and amplified via PCR reaction for 15 cycles. Amplified products underwent final purification using 0.8X Agencourt AMPure XP beads and were resuspended in 32.5 µl of resuspension buffer; 30 µl of purified library were recovered for each sample. Libraries were quantified on the Agilent BioMek Synergy H1 Plate reader with Picogreen reagent. Quality of completed libraries was assayed on Agilent 4200Tapestation with HS DNA1000 ScreenTapes. Libraries were standardized and paired, then 150 bp sequencing was completed using the Illumina NovaSeq6000sequencer.

Quality control of raw sequencing reads was performed using FastQC (version 0.12.1; http://www.bioinformatics.babraham.ac.uk/projects/fastqc/). After, reads were aligned to the human reference genome using the STAR aligner (version 2.7.11a) [15], and gene expression levels were quantified using RSEM (version 1.2.7) [16].

Statistical analysis

All computational analyses were conducted in R (version 4.3.0; R Development Core Team, 2023-04-21). Data was normalized using the Median of Ratios method implemented in DESeq2 (version 1.40.2) [17]. DESeq2 was further employed to identify DE genes between conditions. Significance of differential expression was determined through p-value adjustments using the Benjamini-Hochberg method [18], controlling the False Discovery Rate (FDR) at a threshold of 0.05. For visualization purposes, gene expression was standardized across samples, transforming row values into z-scores. Hierarchical clustering and principal component analysis (PCA) were applied to assess the global patterns of gene expression and to explore the effects of treatment groups, gender, and age.

Enrichment analysis

Upon identifying DE genes, enrichment analysis was conducted using org.Hs.eg.db (version 3.17.0) [19] and clusterProfiler (version 4.8.3) packages [20]. P-values obtained from enrichment tests were adjusted using the Benjamini-Hochberg method, with a significance threshold set at 0.05, to identify the most biologically relevant gene ontology (GO) terms. Further post-hoc enrichment analysis was done utilizing Enrichr software (version 2017b) [21, 22] to investigate KEGG pathway and WikiPathway analyses; these are computed by their respective p-value from the Fisher’s Exact test.

Results

Primary VFF cultures with the designated treatments were grown in vitro. Differentiation of fibroblasts into myofibroblasts was confirmed through immunofluorescence staining for αSMA (Fig. 2). αSMA is a known marker for myofibroblasts as it is a critical component of the stress fibers that develop as a result of fibroblasts becoming activated and differentiating into myofibroblasts [2, 23]. Increased fluorescence was noted in tVFF, followed by moderate and inconsistent fluorescence in sVFF, and minimal fluorescence seen cVFF.

Fig. 2

Expression of myofibroblast marker, alpha-smooth muscle actin (αSMA) across experimental groups. tVFF (a) and sVFF (b) positive for αSMA (green), while cVFF (c) negative for αSMA, all mounted with DAPI (blue). Experiments were performed in triplicates and with negative controls. Photos taken at 20x magnification. Scale bar represents 100 μm

In total, 20,178 coding genes were analyzed from bulk RNAseq. Visual inspection of the PCA plot of all treatment groups (Figure 3a) shows heterogeneity between the cVFF cell lines but homogeneity between nVFF and cVFF. The sVFF group is more similar to cVFF groups, while tVFF groups show more homogeneity and cluster independently from cVFF, sVFF, and nVFF groups. Gene expression profiling considered the top 150 DE genes in one VFF group compared to all other VFF groups (e.g., cVFF versus nVFF, sVFF, and tVFF; Figure 3b).

Fig. 3

Gene expression comparisons across all treatment groups. A Principal component analysis (PCA) plot for all treatment groups, demonstrating similarities between cVFF and corresponding nVFF, and sVFF clustering more closely with cVFF compared to tVFF. There is clustering of tVFF not close to cVFF, nVFF, and sVFF. B Transcriptomic heatmap representing the top 150 differentially expressed genes in each experimental group identified by RNA sequencing, where a corresponds with sVFF, b corresponds with cVFF and nVFF, and c corresponds with tVFF. Adjusted P < 0.05. Colored hierarchical clusters across the top panel represent closely related genes based on sequencing data. Red cluster represents genes most upregulated in tVFF, green cluster represents genes most upregulated in cVFF and nVFF, and blue cluster represents genes most upregulated in sVFF

Enriched GO terms were analyzed using the top 150 DE genes in each group (Table 1). No DE genes were identified between cVFF and nVFF, indicating that differences in the tVFF group are due to TGF-β1 treatment. Additionally, neutralization of the sVFF line did not result in significant differences when compared to the untreated sVFF line.

Table 1 Gene ontology (GO) enrichment terms for the top 150 differentially expressed (DE) genes in each experimental group identified by RNA sequencing

Comparisons between cVFF and sVFF yielded 67 DE genes; 44 were downregulated and 23 were upregulated in sVFF (Fig. 4). Genes related to signaling receptor regulator activity, receptor ligand activity, and signaling receptor activator activity were differentially expressed, including CCL2, DPP4, FGF5, IGF2BP1, THNSL2, and WNT16. There were no significant biological pathways related to the up- or downregulated DE genes in the sVFF group when analyzed separately.

Fig. 4

Differential gene expression for sVFF versus cVFF identified by RNA sequencing. Transcriptomic volcano plot exhibiting 67 differentially expressed genes, represented by red points on the plot. 44 were downregulated (lower log₂ fold change) and 23 were upregulated (higher log₂ fold change) in sVFF (adjusted P < 0.05)

Comparisons between tVFF and cVFF yielded 4772 DE genes; of which, 2761 are downregulated and 2011 are upregulated in the tVFF (Figure 5a). DE genes between these groups were further subdivided based on log₂ fold change extremes, including genes with log2 fold change > 1.5 or <-1.5 (Figure 5b).

Fig. 5

Differential gene expression for tVFF versus cVFF identified by RNA sequencing. A Transcriptomic volcano plot exhibiting 4772 differentially expressed genes, represented by red points on the plot. 2761 were downregulated (lower log₂ fold change) and 2011 were upregulated (higher log₂ fold change) in tVFF (adjusted P < 0.05) B Subset of differential gene expression for cVFF (top 4 rows) versus tVFF (bottom 4 rows) including genes with log₂ fold change > 1.5 or <-1.5. Transcriptomic heatmap exhibiting clustering of differentially expressed genes (adjusted P < 0.05). Colored hierarchical clusters across the top panel represent closely related genes based on sequencing data. The first five colored clusters represent genes most upregulated in cVFF, and the last four colored clusters represent genes most upregulated in tVFF

When considering only the top 50 DE genes, nine pathways related to voltage-gated potassium and ion channels were enriched. Genes related to the fibronectin binding pathway, including TNFAIP6, IGFBP5, and ITGB8, are noted in the top 50 DE genes between these two groups. When all DE genes were considered, there are substantially more enriched GO terms. Given the large number of DE genes, DE genes underwent cluster analysis and were grouped into nine distinct clusters, where the first five clusters highlight genes more upregulated in cVFF and the remaining four clusters highlight genes more upregulated in tVFF. Significant GO terms within each cluster were then identified (Table 2). Of note, these clusters highlight individual cell line differences and thus significant GO terms may represent that of the single cell line more than the VFF group.

Table 2 Gene ontology (GO) enrichment terms for differentially expressed (DE) genes between tVFF and cVFF identified by RNA sequencing

Comparisons between sVFF and tVFF yielded 4876 DE genes; of which, 2351 are downregulated and 2525 are upregulated in sVFF. Given the use of tVFF as a scar-like model, the most highly co-expressed genes were of most interest. 94 genes were identified as significantly co-expressed in sVFF and tVFF compared to cVFF (Fig. 6). Within these 94 co-expressed genes, the only significant GO terms were specifically related to the binding of actin cytoskeleton: actin filament binding and actin binding. These enriched GO terms have significant gene overlap, and are driven by genes FLNB, MYO10, and CNN and TPM isoforms.

Fig. 6

Co-upregulated gene expression for sVFF (top row) and tVFF (following 4 rows) versus cVFF (bottom 4 rows) identified by RNA sequencing. Transcriptomic heatmap exhibiting clustering of 94 significantly upregulated (red) genes in both sVFF and tVFF compared to cVFF (Adjusted P < 0.05)

With the extensive DE genes found between tVFF and sVFF, significant GO terms were analyzed (Table 3). Interestingly, significant GO terms related to ECM structure and cell migration were found in upregulated genes in sVFF but not in tVFF. Of note, despite the direct TGF-β1 treatment in tVFF, the TGF-β receptor signaling pathway is enriched by those upregulated genes in the sVFF group, but not the tVFF group. DE genes that were downregulated in the sVFF group (and subsequently upregulated in the tVFF group) resulted in significant GO terms related to transcription and translation.

Table 3 Gene ontology (GO) and signaling enrichment terms for the differentially expressed (DE) genes between tVFF and sVFF identified by RNA sequencing

Discussion

Bulk RNAseq resulted in DE genes across all groups of VFF other than between nVFF and cVFF. When comparing all VFF groups, GO terms reveal that DE genes in tVFF were enriched in collagen binding, integrin binding, and growth factor activity compared to both sVFF and cVFF (Table 1). Enrichment of these biological processes in VFF treated with TGF-β1 is unsurprising given that myofibroblast differentiation has been shown to stimulate collagen production [1] and corroborates reported enrichment found between untreated VFF and TGFβ1-induced myofibroblasts by Friedman et al. (2025) [10]. With increased collagen deposition from myofibroblasts, increased binding of integrins, cell surface receptors that bind to specific extracellular ligand(s) such as collagen or other ECM components, can occur. Of note, integrin binding can control physiological cellular functions, including differentiation, through cytoskeletal organization [24]. Thus, further myofibroblast differentiation can be promoted by integrins mediating the indirect interaction of ECM and intracellular actin cytoskeleton. Furthermore, with TGF-β1 treatment, it is expected that growth factor activity would be increased. Beyond the introduction of the growth factor to in vitro conditions, there are several possible explanations for increased growth factor activity related to the differentiation into myofibroblasts and the progression of wound healing. For example, with acute inflammation and injury, immune cells are stimulated to release cytokines and growth factors, in turn causing surrounding cells to do the same. Evidence also suggests that signaling between integrins, growth factor receptors, and other proteins may inform the progression of fibrosis [25].

When analyzing gene expression between two VFF conditions, comparisons between nVFF and cVFF showed no DE genes. This lack of significance indicates that the cell culture protocol alone, including the use of a low serum media, did not result in any significant change in gene expression, even with the addition of a TGF-β blocking antibody. Additionally, even when the sVFF were treated with the blocking antibody, there were no significant changes to gene expression in these cells. Therefore, the changes in gene expression seen in the tVFF group – which were the control VFF treated with TGF-β1 for five days in culture – can be attributed to the TGF-β1 treatment, and not additional influence of in vitro conditions.

The number of DE genes between sVFF and cVFF was low compared to other groups’ comparisons. It is possible that because the cells were analyzed at passage 5, there was a convergence of gene expression that might not be seen at an earlier passage. Similarity of scarred VFF and normal VFF phenotypes and gene expression at passages beyond passage 1 has been previously observed and reported [26]. However, TGF-β1 treatment of fibroblasts has previously been reported to primarily induce differentiation into myofibroblasts through enrichment in pathways related to actin filament assembly [27]; the current study is the first known study to look at TGF-β1-treated fibroblasts compared to pathological scarred fibroblasts. DE genes that were found may be indicative of more potent differences that can be seen between sVFF and cVFF. sVFF demonstrated DE genes related to signaling receptor activity and receptor ligand activity, which may facilitate the transcriptional cascade leading to scarring pathophysiology.

As mentioned previously, cytokines, a subset of growth factors, are released by activated immune cells – these cytokines are usually pro-inflammatory, resulting in further inflammation and cytokine production. These signals can induce myofibroblast proliferation and binding to ECM ligands to facilitate remodeling, which, in excess, can ultimately lead to scarring [28]. Proteoglycans, such as heparin, contain GAG chains and are components of the ECM scaffold that are produced by myofibroblasts. ECM components normally sequester growth factors and cytokines, but these reserves are carefully released upon injury [29]. Enrichment of ECM- and growth factor-related pathways in sVFF may indicate the occurrence of the inflammatory phase of wound healing within this sVFF cell line.

There were significantly more DE genes found between tVFF and cVFF. Analysis of the top 50 DE genes resulted in genes predominantly found in the fibronectin binding pathway. Fibronectin is another ECM protein that myofibroblasts produce and assemble into a matrix for wound healing faster than inactivated fibroblasts [30]. Fibrillogenesis of collagen requires fibronectin assembly, thus myofibroblasts likely accelerate fibronectin synthesis in order to increase ECM formation during scarring. Additionally, the top 50 DE genes were enriched in nine pathways related to channel activity. Recent literature has documented the role of ion channels in myofibroblast differentiation and fibrosis throughout the body, including the kidneys [31], lungs [32], and heart [33, 34]. More specifically, TGF-β1 increases susceptibility to atrial fibrillation, and this is attributed to the effect of TGF-β1 on mRNA expression and subsequent encoding of calcium channels [35]. A model of crosstalk between cation channels and the TGF-β signaling pathway has also been described in pulmonary fibrosis [32], which has led to a focus on ion channels as potential therapeutic targets in fibrosis [36,37,38].

Significant GO terms resulting from upregulated genes in tVFF include growth factor activity, growth factor receptor binding, and transmembrane receptor protein tyrosine kinase activator activity, among others (Table 2). Given that tVFF were treated with TGF-β1, this may play a significant role in the enrichment related to growth factor activity and binding; however, this treatment may also inform the enrichment of the transmembrane receptor protein tyrosine kinase activator activity. TGF-β1 may bind to any of the three isoforms of the TGF-β receptor (TGFBR), and two of the three are tyrosine kinases [39]. These kinases are necessary for activating downstream pathways, including the TGF-β signaling pathway. A similar enrichment pattern can be appreciated within sVFF, where signaling receptor activity was also enriched. This enrichment may also be related to downstream signaling from TGF-β or similar, despite not receiving the direct TGF-β1 treatment.

Finally, while there was a very large number of DE genes between tVFF and sVFF, a crucial similarity was seen: enrichment of the actin binding pathway in both conditions. This enrichment is a notable indication of activation in TGFβ1-induced myofibroblasts [10]. Actin stress fibers are a defining morphologic feature of myofibroblasts, playing a critical role in the contractile function necessary for these activated cells to close a wound [40]. αSMA is incorporated into these stress fibers to add contractile force, and focal adhesions (another enriched process found between these groups) serve as a “home-base” for the formation of these actin fibers [40]. Despite the large number of DE genes, analysis of co-upregulation reveals that both tVFF and sVFF achieve myofibroblast differentiation and pathophysiology. This further supports the use of tVFF as a model of a “scar-like” phenotype in in vitro experiments. Because only one line of sVFF was available in this investigation, comparisons between tVFF and sVFF were focused mainly on the indicators of myofibroblast differentiation rather than investigating novel differences between the experimental groups.

Several limitations of this work are worth noting. One caveat is that our results are strictly based on VFF transcription, in vitro. Despite care to maintain these cells at very low passages, gene expression from in vivo fibroblasts may display differences due to interactions with other cells and systemic influence. Additionally, as previously highlighted, VFF have both unique positioning and exposure to mechanical stress and environmental irritants [5]. While this may be advantageous in order to study transcriptional response and resulting pathology from such exposure, gene expression of fibroblasts in areas of the body that have different position and exposure (e.g., liver or kidney) may exhibit significantly different transcriptomics, both normally and pathologically. VFF may better model transcriptomics of fibroblasts found in the lung or heart due to more similarities in biomechanics. Finally, access to fibroblasts derived from vocal fold pathologies is significantly limited due to the difficulty in retrieving these cells without causing further damage to the vocal fold. The use of the sVFF line should be considered as a comparator to tVFF to corroborate the differentiation into myofibroblasts; transcriptional signatures (e.g., enrichment of actin binding and assembly) seen in both the sVFF line and tVFF provide an indication that myofibroblast differentiation occurred.

Conclusions

Fibroblasts remodel and maintain the ECM in homeostasis; [41] upon differentiation into myofibroblasts, ECM production increases to aid in contraction and cell signaling changes to favor wound healing. As a result, transcriptomics of normal VFF understandably differ from myofibroblasts, including from those retrieved from a vocal fold scar and those treated with TGF-β1 to model myofibroblasts from scarring. This is the first research to investigate these transcriptional differences, in bulk, between normal and pathological fibroblasts and further understand how these cells change and respond in injury. Upregulated genes in cVFF were enriched in pathways related to cell adhesion and cell cycle/division, while upregulated genes in tVFF and sVFF were enriched in pathways related to ECM binding and growth factor activity. Both tVFF and sVFF demonstrate enrichment in actin binding, a defining difference between fibroblasts and myofibroblasts. Despite large differences in transcriptomics between tVFF and sVFF, tVFF serve as a useful in vitro model of myofibroblasts and highlight key similarities to myofibroblasts extracted from scar pathology, as well as expected differences related to normal fibroblasts from healthy vocal folds.