RNA sequencing identifies transcriptional changes in the rabbit larynx in response to low humidity challenge

Background Voice disorders are a worldwide problem impacting human health, particularly for occupational voice users. Avoidance of surface dehydration is commonly prescribed as a protective factor against the development of dysphonia. The available literature inconclusively supports this practice and a biological mechanism for how surface dehydration of the laryngeal tissue affects voice has not been described. In this study, we used an in vivo male New Zealand white rabbit model to elucidate biological changes based on gene expression within the vocal folds from surface dehydration. Surface dehydration was induced by exposure to low humidity air (18.6% + 4.3%) for 8 h. Exposure to moderate humidity (43.0% + 4.3%) served as the control condition. Ilumina-based RNA sequencing was performed and used for transcriptome analysis with validation by RT-qPCR. Results There were 103 statistically significant differentially expressed genes identified through Cuffdiff with 61 genes meeting significance by both false discovery rate and fold change. Functional annotation enrichment and predicted protein interaction mapping showed enrichment of various loci, including cellular stress and inflammatory response, ciliary function, and keratinocyte development. Eight genes were selected for RT-qPCR validation. Matrix metalloproteinase 12 (MMP12) and macrophage cationic peptide 1 (MCP1) were significantly upregulated and an epithelial chloride channel protein (ECCP) was significantly downregulated after surface dehydration by RNA-Seq and RT-qPCR. Suprabasin (SPBN) and zinc activated cationic channel (ZACN) were marginally, but non-significantly down- and upregulated as evidenced by RT-qPCR, respectively. Conclusions The data together support the notion that surface dehydration induces physiological changes in the vocal folds and justifies targeted analysis to further explore the underlying biology of compensatory fluid/ion flux and inflammatory mediators in response to airway surface dehydration.


Background
Voice disorders are a prevalent communication disorder affecting human health worldwide (1)(2)(3)(4)(5)(6). In the United States general population, the prevalence of voice disorders has been estimated at 6.2% (7), and more recently, at 7.6% (8). Data from the National Longitudinal Study of Adolescent to Adult Health shows the same 6% estimate among the adolescent population (9). The development of voice disorders is identi ed as an occupational hazard, particularly among speakers who depend on a healthy voice for their livelihood. School teachers, entertainers, legal professionals are all at greater risk of dysphonia from voice disorders (3, 7, 10-13). The economic impact of voice disorders is substantial. The average associated health care costs in the United States have been estimated at almost 200 million dollars (14), and a study of Brazilian teachers having to take time away from work due to dysphonia illustrates the potential impact of a loss of productivity in the workforce (15). Taken together, the impact of voice disorders on society supports the need for a more comprehensive understanding of the development of voice disorders and therapies to address them.
Interventions for voice disorders exist along a continuum of non-invasive behavioral modi cations to phonosurgery. The focus of this study is restricted to the consideration of avoiding dehydration of the laryngeal surface as a commonly prescribed prophylactic and therapeutic practice among speechlanguage pathologists (4,(16)(17)(18)(19). An absence of clear mechanistic physiological associations between dehydration and dysphonia challenges the validity of currently prescribed treatment protocols and precludes informing the development of new therapies.
Dehydration, as it relates to voice, occurs under two paradigms: systemic dehydration and airway surface dehydration. Systemic dehydration, decreased total body water, has been shown to negatively impact phonatory effort in humans and acoustic measures in humans and ex vivo animal models (20)(21)(22)(23). The presumed mechanism described in the literature is changes to the biomechanical properties of the vocal folds due to decreased water (4,(24)(25)(26)(27)(28). Surface dehydration as related to voice is de ned as loss of water from the luminal surface of the larynx and vocal folds. In everyday life, this may be caused by exposure to air of low humidity or increased respiratory rate from exercise. While there is evidence suggesting that surface dehydration within the larynx negatively impacts phonation with similar outcomes as systemic dehydration, recent studies in humans (29)(30)(31)(32) do not always nd a signi cant correlation between the two. Further, the functional nature of many investigations leaves many questions about the underlying biology unanswered.
Unfortunately, rigorous in vivo analysis of the physiology of laryngeal surface dehydration is precluded by the invasive nature of data collection and the ethical implications of causing vocal injury in human subjects. Human studies are, therefore, generally limited to acoustic and aerodynamic measures or postmortem evaluation. Conversely, animal models have largely allowed for ex vivo studies, which provide ample evidence that surface dehydration impacts vocal fold biomechanics and function (33)(34)(35), but the molecular pathobiology and resulting homeostatic compensatory mechanisms remain unclear. An attractive surrogate to the vocal folds is the airway distal to the larynx, which has been studied in the context of airway surface uid homeostasis and response to luminal perturbations (36)(37)(38). It has long been established that the humidity of inspired air can affect the magnitude of water lost to respiration (39) and that the resulting concentration of luminal electrolytes can cause dramatic physiological responses in the trachea, upper and lower airways (40,41). The vocal folds are covered by nonkeratinized strati ed squamous epithelium, and the laryngeal lumen is predominately covered by respiratory epithelium. Therefore, the larynx may respond to perturbations similarly to the tracheal epithelium. This potential is supported in studies assessing vocal fold ion ux to altered composition of luminal surface uid (42)(43)(44). However, these were in vitro studies limiting the generalization of the data. Further studies are required to address questions of the speci c underlying biology.
To probe for potential physiological responses to surface dehydration, we used an in vivo rabbit model. Anatomically, the rabbit larynx is grossly similar to the human larynx. Its size has been approximated to 8.6 × 5.5 mm at the level of the arytenoids (45,46), consistent with the dimensions of the human newborn larynx (47). Additionally, the literature demonstrates that rabbit larynges exhibit su cient biological similarity to humans and have been used in molecular and histological studies of the vocal folds (48)(49)(50)(51)(52). The rabbit larynx has also been used to characterize the physiological response to injury secondary to phonation (51,52) or laryngeal and vocal fold surgery (53)(54)(55). The common use of rabbits for laryngeal studies and the relative small size for handling and housing makes this animal a suitable model for this study.
In this study, we sought to identify transcriptional-level changes in response to low humidity exposure that suggest a physiological response to surface dehydration within the membranous vocal folds or the vocal fold lamina propria. We successfully addressed the following aims: [1] construction and evaluation of an environmental chamber capable of exposing rabbits to a consistent, physiologically-realistic low relative humidity environment and [2] investigation of the effects of 8 hours of low humidity exposure on rabbit larynx by way of RNA sequencing (RNA-Seq). An 8-hour exposure was selected as representative of a typical working day for human subjects. We used low humidity rather than desiccated air as the surface dehydration challenge to increase the ecological validity of the study. Rabbits exposed to moderate humidity served as the control condition.

Humidity Challenge and Gross Physical Assessment
A total of eight rabbits were challenged with low humidity, and six rabbits were exposed to moderate humidity (control condition) in a specially fabricated environmental chamber ( Fig. 1; see Methods section for details). Low humidity was 18.6 ± 4.3% (mean ± standard deviation) over the 8 hours. The moderate humidity exposure was 43.0 ± 4.3% over the 8 hours (Fig. 2). There was no observable behavioral differences or evidence of respiratory distress following exposure in either group. No gross evidence of in ammation or damage to the laryngeal mucosa was observed during visual examination under a dissecting microscope.

Packed Cell Volume (PCV)
The pre-experiment PCV (%) across all 14 rabbits was 46.7 ± 2.8 (mean ± SD). The % change in PCV from baseline to after the experiment did not differ signi cantly between the low and moderate humidity groups (p = 0.1692).

Sequence Read Mapping and RNA-Seq
Approximately 69 to 112 million paired reads were obtained by RNA-Seq with an average of 70% quality reads mapping to genes in the rabbit genome. Differential gene expression by Cuffdiff revealed 103 genes reaching an FDR ≤ 0.5 with 61 meeting the additional fold change (|log2 FC| ≥ 1) ltering criterion. Of these, 48 genes were considered signi cantly downregulated and 13 genes were signi cantly upregulated. The 10 genes with the greatest up-and downregulated fold changes from this list of 61 are shown in Table 1. A complete list of all genes identi ed is provided within the Additional le 1: Table S1.

Functional Enrichment Analysis
Functional enrichment analysis by DAVID and STRING provided similar but distinct sets with FDR < 0.05.
DAVID identi ed 5 GO terms for biological process, 7 GO terms for cellular component, 4 GO terms for molecular function, and 7 processes by KEGG. Genes clustered into three groups, including KEGG processes related to cardiac muscle function, chemical carcinogenesis, and ECM-receptor interaction. STRING provided a richer set with 7, 15, and 19 GO terms for biological process, cellular component, and molecular function, respectively, and 2 KEGG processes. The full lists of terms, functions, associated genes, and statistics for the aforementioned analyses are provided in Additional le 2: Table S2.
The predicted protein-interacting network generated by STRING is shown in Fig. 3. There were 8 clusters identi ed with between 2 to 10 gene products. Larger clusters contain members that are associated with cellular response to external stimuli and immune response (dark green, lavender), muscle function (red), keratinocyte development (light green), and ciliary function (aqua).

Discussion
The transcriptional changes observed in this study indicate that just 8 hour exposure to a low humidity environment can adversely affect vocal fold biology. To the best of our knowledge, this is the rst study to demonstrate the effects of surface dehydration on vocal fold tissue in vivo. Surface dehydration was induced in a physiologically-realistic manner by exposure to low humidity. It is noteworthy that exposure to low relative humidity below the Occupational Health and Safety Administration (OHSA) recommended limit of 20% induced transcriptional changes within functional gene categories including in ammation, ion transport, and keratinocyte development.
The most robust functional enrichments identi ed by STRING were stress, defense, and in ammatory responses that included downregulated ORM1, S100A9, S100A12, and SAA1, and upregulated IL1RN, MCP1/2, and MMP12 genes. Additionally, outside of the STRING analysis, various genes for immunoglobulin chains were identi ed, three of which were downregulated and one that was upregulated. Interestingly, this cluster presents two opposing interpretations of innate immune dampening and possible macrophage activation.
While none of these genes or corresponding proteins are described within the larynx, the downregulated cluster can be interpreted as a dampening of acute in ammatory response. ORM1 and SAA1 are both acute phase proteins. ORM1 is an acute phase protein that has been shown to polarize M2 macrophage differentiation (63) and to enhance epithelial integrity in a culture model of the blood-brain barrier (64). While ORM1 exhibits anti-in ammatory activity and its downregulation may allow for the development of a more robust in ammatory process, it may also be interpreted as indicative of surface dehydration not contributing to an activating in ammatory event. SAA1 is also an acute phase protein and is associated with a variety of pathological conditions, but it has also been shown to positively in uence keratinocyte activity (65). The S100 proteins are diverse with involvement in several cellular processes, but both S100A9 (66) and S100A12 (67) have been described as damage associated molecular patterns in the literature. Taken together, these results suggest that either surface dehydration is not inducing in ammatory pathways or that there is active repression of pro-in ammatory mediators. The latter is substantiated by the increase of IL1RN which encodes the IL-1 receptor antagonist (IL1RA). IL1RN was upregulated in the posterior cricoarytenoid muscle one week following transection of the recurrently laryngeal nerve in a rat model (68), and IL1RA was signi cantly increased following 8 hours of industrial exposure to respirable and inhalable dust in humans (69). Together this substantiates a role for the increased IL1RN we observed and of a possible active innate immunity repression in response to the low humidity challenge.
Conversely, the upregulation of MMP12 and MCP1 genes may suggest the activation of in ammatory macrophages. MMP12 was the most signi cantly upregulated gene in this study by RNA-Seq and RT-qPCR. MMP12 exhibits proteolytic activity on multiple ECM components including elastin, bronectin, entactin, and type IV collagen (70), all of which are expressed within the vocal folds. Although called "macrophage elastase", it is also expressed in human vocal fold broblasts (71) and bronchial epithelial cells in vitro (72), and in both super cial and deep epidermal layers of the skin in response to ultraviolet radiation (73). MMP12 has a potential role in the development of dysphonia following low humidity exposure since type IV collagen and elastin play an important role in the viscoelasticity and phonatory function of the vocal folds (74,75). MMP12 may contribute directly to in ammation though epidermal growth factor receptor (EGFR) dependent induction of IL-8 from the respiratory epithelium (76). Interestingly, MMP12 has been shown to positively in uence wound healing following epithelial injury to the cornea (77), so it is unclear if the upregulated response to low humidity would be deleterious or in uence a reparative response in the vocal folds. MCP1 is an α-defensin expressed in the lungs of fetal and adult rabbits (78); it is secreted from neutrophils and rabbit lung macrophages and exhibits broad antimicrobial activity In our study, the expression of MCP1 was novelly detected in the rabbit larynx, and its upregulation in repsose to low humidity warrants further investiation including targeted anaylsis of differential expression between in ammatory cells and the larygeal tissue.
It is not surprising to nd evidence of a pro-in ammatory response with surface dehydration as other environmental stressors such as simulated acidic re ux (80), hypertonic challenge (43), and phonotrauma (52,81) can perturb the epithelial tight junctions of the vocal folds-indicative of the activation of proin ammatory pathways. As we did not investigate for cell-speci c gene expression in this study, we are limited to conclude if the upregulation of these genes re ects activation of macrophages or activity of the epithelium or lamina propria broblasts, and further study is warranted.
An intriguing hypothesis for a case of macrophage activation would be altered response to local microbiome or pathogens resulting from changes to the laryngeal microenvironment following dehydration.
The perturbation of ion transport or other lubrication mechanisms is anticipated as a response to the altered hydration state of the laryngeal surface (82). Although no gene or protein interaction enrichment cluster was identi ed within the 103 DEGs analyzed, presumably due to the diversity of substrate and transporter type, a considerable set of ion and solute transporter related genes were identi ed by RNA-Seq, including ECCP, SLC5A1, SLC13A5, SLC23A1, SLC27A2, and ZACN. All SLC family members were downregulated. This set represents predominantly ion transport, with SLC13A5 and SLC27A2 being involved in glucose transport and fatty acid ligation. In vitro studies of human nasal epithelial cells (41) and human bronchial cell culture (83) demonstrated that apical osmotic pressure can result in altered epithelial electrolyte transport; however, studies with canine tracheal and bronchial cell culture (84) and an in vivo canine model (85) concluded that not all epithelial uid ux is coupled to electrolyte transport. This evidence suggests that the epithelium may respond to either aberrant electrolyte concentrations or non-ionic osmotic pressure. It is not surprising to nd evidence of altered chloride secretion speci cally, as balanced sodium and chloride ion secretion is attributed to volume regulation of the airway surface uid, but the contribution of transport of other ionic and non-ionic species is not well described for airway surface uid regulation. Our results suggest the pertinence of future targeted study of noncanonical secretion products in the respiratory tract.
Although the ECCP is annotated as an epithelial chloride channel protein, the translation product for ECCP is neither well characterized nor has a direct ortholog in humans. It may belong to the calciumactivated chloride channel proteins (CLCA) family as identi ed by conserved functional domains, although it exhibits limited homology to the rabbit CLCA proteins. The genes for CLCA1, CLCA2, and CLCA4 lie within the same genetic neighborhood as ECCP but were identi ed by RNA-Seq with FDR > 0.99, indicating they are not differentially expressed in our model of surface dehydration (Additional le 1: Table S1). This suggests a distinct role for ECCP and its downregulation that warrant further investigation as an ion channel protein newly described in the context laryngeal surface dehydration. In contrast to ECCP, ZACN was upregulated in low humidity compared to moderate humidity but failed to reach statistical signi cance by RT-qPCR. ZACN is a cation channel expressed in the human trachea and other tissues and demonstrates permeability to potassium ions but not to chloride ions (86); there is no discussion of its expression in the vocal folds in the literature, and it is unclear if it may also be sodium ion permeable. Taken together with the SLC family members identi ed, these results support a potential role for solute ux as a homeostatic response to surface dehydration. Interestingly, however, the downregulation of chloride transportation would be a counterintuitive response to surface dehydration at the apical membrane as chloride is generally directed out of the cell and aberrant chloride transport can be detrimental in the airways as seen in cystic brosis. There is a distinction between the respiratory epithelium of the airways and the nonkeratinized strati ed squamous epithelium of the vocal folds, so care must be taken with direct translations of actions between the two.
The mucins are equally important to maintain satisfactory hydration of the laryngeal surface as ion and uid ux. MUC12, MUC21, and TFF1 were identi ed as downregulated by RNA-SEq. Both mucins are members of the cell-surface associated mucin family, and as such, should originate directly from the epithelial cells. The rst exon of MUC12 exhibited increased expression in laryngeal epithelium from laryngeal re ux patients compared to re ux negative patients (87). Exogenous surface expression of MUC21 in in vitro cell culture reduced intercellular adhesion and adhesion to extracellular components (88). It is interesting then to observe all three to be downregulated. However, in addition to roles as epithelial protectants, mucins and related proteins also serve roles in cell signaling with physiological consequences. This is recently shown for MUC21 overexpression as in uencing the development of lung adenocarcinoma (89) and TFF1 in uencing epithelial-mesenchymal transition. Together, this may be a contributing factor to the STRING cluster of keratinocyte differentiation factors discussed below, but further study is warranted to determine which cell types are expressing these genes and which cell signaling may be impacted.
Although there was no gross in ammation observed, some level of epithelial cellular response to surface dehydration is expected. The vocal folds are covered by a non-keratinized strati ed squamous epithelium for which some aspects of development are well understood, such as embryological developmental factors and differentially expressed structural components (90, 91), but a comprehensive molecular description is not available as for other epithelia like the epidermis. It is interesting that several keratinocyte developmental factors were identi ed with RNA-Seq and as a protein interaction cluster in the STRING analysis: CDSN, CNFN, CRNN, KRT80, KRTDAP, and TGM3. Also identi ed by RNA-Seq were SPBN, another keratinocyte factor, and CDHR4, a cell interaction mediator. All of these were downregulated. SPBN, CDHR4, and CDSN were selected for RT-qPCR validation. All three gene products may be involved in maintaining the integrity of the strati ed squamous epithelium, though none have been described speci cally within the vocal folds until this study. SPBN is expressed in the suprabasal layers of tongue, stomach, and epidermis (92). It is required for keratinocyte differentiation in an in vitro skin model (93) and skin development in murine embryos (94). The speci c activity of CDHR4 is not described in the literature, but family member CDHR2 is expressed in gastrointestinal epithelial cells and is associated with microvillus development (95), while family member CHDR3 is expressed in ciliated respiratory epithelial cells and is associated with ciliary development and intercellular interactions (96). CDSN is expressed in the stratum granulosum of human skin and appears to participate in cellular cohesion at this level, with its loss associated with desquamation (97,98). That the entire cluster was downregulated substantiates surface dehydration as capable to in uence vocal fold epithelial maintenance. Further study is required to elucidate the speci c roles of these proteins within the vocal folds, as this epithelium is distinct from the epidermis.
As part of this study, we developed a method to e ciently challenge rabbits to low humidity.

Limitations
A limitation of designing an environmental chamber as described here was that it precluded the provision of relative humidity lower than 15%. While environmental rooms and chambers are commercially available, they are cost prohibitive and their small size precludes the use of certain animal models, such as rabbits. Another limitation of the study is that we only observed a single time point after low humidity exposure. It has been shown that local response to vocal fold injury is transient and time-dependent (48,54,100). Further studies speci cally observing for in ammatory response at multiple times points within a single challenge or within repeated or chronic challenges would be helpful in further characterization of vocal fold biology. Finally, the dissected vocal fold tissue included striated muscle and small amounts of respiratory epithelium immediately above and below the region of the vocal folds. Therefore, genes associated with muscle or respiratory epithelium were not selected for the discussion.

Conclusions
In this study, we investigated whether surface dehydration induced by low humidity would affect vocal fold biology. We successfully developed an e cient and cost-effective environmental chamber to induce surface dehydration. Humidity was maintained at occupationally relevant levels to observe for physiological impact. Low humidity challenge did not induce systemic dehydration as evaluated with PCV. High-throughput RNA-Seq provided evidence of a biologic response to low humidity challenge, including changes within stimulus-response, lubricative mechanisms, and potential alterations in epithelial maintenance. This study successfully serves as a framework for targeted investigations of molecular response to surface dehydration in the larynx.

Animals
All experiments were conducted in accordance with the guidelines and after approval of the Purdue Animal Care and Use Committee (Protocol # 1606001428). Animals were obtained from Envigo Global (Indianapolis, IN) and acclimatized for at least one week. Male New Zealand White rabbits, six to eight months of age, and approximately 3 Kg were used for all experiments. Each experiment was run with two rabbits at a time randomly assigned to either the low (n = 8) or moderate (n = 6) humidity group. Samples sizes were selected based on recommendation from the Purdue Bioinformatics Core to ensure ideal minimum samples for statistical validity of RNA-Seq (n = 6 from each group). Changes to PCV were examined for all rabbits. RT-qPCR validation was conducted with 13 rabbits (low n = 7, moderate n = 6); one rabbits was excluded due to poor quality of RNA obtained after repeat extraction. Food and water were withheld during experiments under both humidity conditions. To encourage consistent, baseline hydration, all animals were pre-hydrated with 0.1 M sucrose in water ad libitum for the two days preceding the experiment. Euthanasia was completed by a single 1.0 mL IV dose of Beuthanasia-D Special (Schering Plough Animal Health Corp., Union, NJ).

Humidity Challenge Protocol
Eight hour low humidity and moderate humidity exposure were conducted in a specially fabricated environmental chamber. The chamber interior was segmented into three similar compartments, each with dimensions approximately 61 cm × 61 cm × 46 cm (Fig. 1a, b). Two compartments were sealed to limit the in ux of room air and were intended for low humidity exposure, whereas the third compartment was left open to room air and was intended for a moderate humidity control. Gated duct caps were included within the wall of the low humidity compartment to allow for titration of room air as necessary.
Low humidity was achieved with a 70-pint commercial dehumidi er (Hisense DH70K1G: Qingdao, China) set to High Continuous attached to the chamber via 4-inch ducting. Moderate humidity exposure was achieved by opening the chamber airspace to room air without conditioning from the dehumidi er.
Internal relative humidity and temperature were tracked using a HOBO Data Logger with a 12-bit Temperature/Relative Humidity Smart Sensor (U14-002, S-THB-M002: ONSET, Bourne, MA) at one-minute intervals.

Blood Collection and Analysis
Blood was collected in heparinized tubes at the beginning of the 8-hour experiment and immediately prior to euthanasia via venipuncture of the lateral ear vein to minimize trauma and distress of collection.
Packed cell volume (PCV) was measured manually by visual assessment using a microhematocrit reader card following centrifugation.

Sample Collection and RNA Extraction
The larynx and proximal trachea were excised from each animal immediately following euthanasia. The larynx was bisected posteriorly along the sagittal midline and pinned onto wax to expose the laryngeal lumen. Counts from all replicates were merged using custom Perl scripts to generate a read count matrix for all samples.

Differential Gene Expression Analysis and Annotation
Differential gene expression analysis between low and moderate-humidity groups was carried out using 'R' (v 3.

Functional Enrichment Analysis and Predicted Protein Interactions
Gene expression data were ltered by two levels: a false discovery rate (FDR) of less than or equal to 0.05, and a log-2 fold change (FC) of expression with an absolute value greater than or equal to 1.0. FC positive values imply upregulation, while negative values imply downregulation of genes in the vocal folds exposed to low humidity versus moderate humidity challenge. The set of genes meeting both criteria were used for functional enrichment analysis while the full set of genes meeting the FDR criterion was used for prediction protein interaction analysis.
GO and KEGG enrichment analyses were performed using both DAVID (v6.8) (59) and STRING (Search Tool for the Retrieval of Interacting Genes/Proteins) v.11.0 (https://string-db.org). STRING was also used to generate a network of predicted protein interactions based on the same data. STRING analysis parameters were set with line thickness representing the strength of the data to support interaction including text mining, experimental, database, co-expression, neighborhood, gene fusion, and cooccurrence sources, the minimum required interaction score set to 0.4, shell parameters set to "None", and disconnected nodes to be hidden. Clusters were generated based on the Markov Cluster Algorithm with the in ation parameter set to 2. All the analyses were performed using the differentially expressed genes Ethics approval and consent to participate All experiments were conducted in accordance with the guidelines of the Purdue Animal Care and Use Committee (Protocol # 1606001428).

Consent for publication
Not applicable.

Availability of data and materials
Sequencing data were submitted to the NCBI GEO database under accession number GSE148588. A full list of genes identi ed and Cu inks-calculated differential expression is available within the supplementary information.

Competing Interests
The authors identity no competing or con icts of interest.

Funding
Funding was provided from R01DC0115545 (National Institutes of Health/National Institute on Deafness and other Communication Disorders). The NIH/NIDCD did not participate in any component of the study design, execution, or analysis.
Author's Contributions TWB was responsible for the conception and design of the study, data acquisition, data analysis and interpretation, and manuscript writing. APS supervised the study, provided resources, reviewed and edited the manuscript. NCN was involved in the design of the study, data analysis, and reviewing and editing the manuscript. MPS was responsible for conception and supervision of the study, funding, reviewing and editing the manuscript. AC: was responsible for conception and supervision of the study, project administration, and reviewing and editing the manuscript. All authors approved the nal manuscript.  Protein interaction network was created using STRING. The 103 differentially expressed genes identi ed by Cuffdiff with an FDR ≤ 0.05 were loaded into the system. The line thickness represents the strength of the data to support the interaction, including text mining, experimental, database, co-expression, neighborhood, gene fusion, and co-occurrence sources. The minimum required interaction score was set to 0.4. Shell parameters were set to "None". Disconnected nodes are not shown. Cluster colors are based on the Markov Cluster Algorithm with the in ation parameter set to 2.