Comprehensive transcriptome analysis of fluid shear stress altered gene expression in renal epithelial cells

Renal epithelial cells are exposed to mechanical forces due to flow‐induced shear stress within the nephrons. Shear stress is altered in renal diseases caused by tubular dilation, obstruction, and hyperfiltration, which occur to compensate for lost nephrons. Fundamental in regulation of shear stress are primary cilia and other mechano‐sensors, and defects in cilia formation and function have profound effects on development and physiology of kidneys and other organs. We applied RNA sequencing to get a comprehensive overview of fluid‐shear regulated genes and pathways in renal epithelial cells. Functional enrichment‐analysis revealed TGF‐β, MAPK, and Wnt signaling as core signaling pathways up‐regulated by shear. Inhibitors of TGF‐β and MAPK/ERK signaling modulate a wide range of mechanosensitive genes, identifying these pathways as master regulators of shear‐induced gene expression. However, the main down‐regulated pathway, that is, JAK/STAT, is independent of TGF‐β and MAPK/ERK. Other up‐regulated cytokine pathways include FGF, HB‐EGF, PDGF, and CXC. Cellular responses to shear are modified at several levels, indicated by altered expression of genes involved in cell‐matrix, cytoskeleton, and glycocalyx remodeling, as well as glycolysis and cholesterol metabolism. Cilia ablation abolished shear induced expression of a subset of genes, but genes involved in TGF‐β, MAPK, and Wnt signaling were hardly affected, suggesting that other mechano‐sensors play a prominent role in the shear stress response of renal epithelial cells. Modulations in signaling due to variations in fluid shear stress are relevant for renal physiology and pathology, as suggested by elevated gene expression at pathological levels of shear stress compared to physiological shear.

It is currently not known in detail how fluid shear stress affects cellular behavior and which signaling pathways are altered. Furthermore, gene expression and the overall cellular behavior will be the effect of an integration of the different signaling pathways, triggered by shear stress and by cytokine stimulation. In this study we set out to obtain a comprehensive overview of the transcriptome under static and shear stress conditions in renal epithelial cells to get more insight in the pathways and processes involved in the shear response.
Therefore, we applied RNA-sequencing as an unbiased means to interrogate renal epithelial cell type-specific transcriptome alterations upon fluid shear stress. Our data indicate that genes involved in TGF-β, MAPK, and Wnt signaling are up-regulated by shear stress, while the JAK-STAT related genes seems to be down-regulated.
Using ALK4/5/7 and MEK1/2 inhibitors, we showed that the shear stress-induced signaling cascades are largely modulated by TGF-β/ ALK5 and MAPK/ERK signaling. Cilia removal abrogated shear induced gene expression of a subset of genes, but genes involved in TGF-β, MAPK, and Wnt signaling were hardly affected, suggesting that other mechano-sensors also play an evident role in the shear stress response of renal epithelial cells. Furthermore, altered expression of genes involved in cell-matrix, cytoskeleton and glycocalyx remodeling, as well as amino acid, carbohydrate, and cholesterol metabolism, indicate that shear stress is regulating gene expression at several levels for cellular homeostasis. Finally, we showed that expression of several genes is elevated at pathological levels of shear stress compared to physiological controls, suggesting that variations in fluid shear stress might be relevant for the pathology in kidney diseases due to an imbalance in cellular signaling. The Netherlands) and ammonium sulfate (#A-2939) from Sigma-Aldrich (Zwijndrecht, The Netherlands) were used as previously described (Kunnen et al., 2017).
Differentially expressed genes were selected with an adjusted p-value (corrected for multiple hypotheses testing) of <0.05. Count per million (CPM) values were calculated by dividing the read counts by total read counts of the sample, which is a measure for the abundance of the transcript. CPM > 2 was used to exclude low expressed genes.

| Quantitative PCR
Gene expression analysis by quantitative PCR (qPCR) was performed as described previously (Happe et al., 2011 were tested using one sample t-tests. One-way analysis of variance (ANOVA) was used when cells were exposed for a different time or to a different flow rate. Two-way analysis of variance (ANOVA) was used, when the shear stress response was compared to a second treatment. The ANOVA was followed by post-hoc Fisher's LSD multiple comparison, if the overall ANOVA F-test was significant.
p < 0.05 was considered to be statistically significant. 3 | RESULTS

| Fluid shear stress induced transcriptional changes in PTECs
To study genome wide fluid-flow induced cellular alterations, proximal tubular epithelial cells (PTEC) were exposed to fluid shear stress of 1.9 dyn/cm 2 using a cone-plate device. Controls were similarly treated under static conditions. After 6 hr fluid shear stress or static exposure, total RNA was isolated and gene expression was analyzed using next generation sequencing (NGS) on the Illumina HiSeq 2500 platform.
After quality checks the reads were aligned to mouse genome (GRCm38) and gene expression was quantified using HTSeq-Count.
Count per million (CPM) values were calculated as a measure for the abundance of the transcript (Supplementary Table S2).  shows increased expression of genes encoding proteins involved in TGF-β ligand activation (Furin, Thbs1) or ligand inhibition (Ltbp2), the transcription factor Smad3, but also the inhibitors Smad7, Smurf1, Skil, and Tgif1, all critical components of the pathway (Table 2).
However, the most prominently activated signaling pathway is the mitogen-activated protein kinase (MAPK) pathway (Tables 2 and   S4 (Tables 3 and S5). This also includes PI3K/AKT related signaling, which is not included as core signaling pathway from KEGG in the MSigDB.  Pathway analysis (KEGG, Reactome and Biocarta) done on 813 significantly up-regulated genes upon fluid shear stress in PTECs using MSigDB. The most significantly altered core signaling pathways are shown and ordered together, followed by the lowest false discovery rate (FDR). K, number of genes in pathway database; k, number of genes in overlap. For the complete list of the pathway analysis of up-regulated genes upon fluid shear stress see Supplementary Table S4.
After 16 hr gene expression was significantly increased for all tested genes (Supplementary Figure S2).
In addition to increased expression of canonical SMAD2/3 targets,

| Primary cilia only play a role in a part of the shear stress response in PTECs
Since defects in cilia formation and function have profound effects on the development and physiology of kidneys and other organs (Goetz & Anderson, 2010;Quinlan et al., 2008), we investigated the shear stress response in PTECs after cilia removal by ammonium sulfate. Expression of Plk2, Prune2, and Ets1 were clearly cilia dependent, since the shear stress induced response was completely lost after cilia ablation (Figure 4). In contrast, genes involved in TGF-β, Wnt, MAPK, and JAK/STAT signaling, that is, Junb, Fbln5, Wisp1, Map3k20, Map4k4 as well as Stat1, were only slightly or not affected in the shear stress response upon cilia removal (Figure 4). Although shear induced down-regulation of Jak2 was abrogated, Jak2 expression in static cells was already reduced upon ammonium sulfate treatment ( Figure 4). Our data suggests that shear stress regulated gene expression in PTECs is only partially cilia dependent and other mechano-sensors are involved as well.

| Shear stress induced gene expression in PTECs is flow rate dependent
Thus far we applied fluid shear stress of 2.0 dyn/cm 2 , which is known to be an increased physio-pathological shear stress (Essig, Terzi, Burtin, & Friedlander, 2001;Grabias & Konstantopoulos, 2012Weinbaum et al., 2010). To compare the gene expression to physiological levels of shear, we exposed the cells to a shear stress range of 0.25-2.0 dyn/cm 2 .
Expression of Wisp1, Map3k20, Map4k4, and Ets1 was clearly flow rate dependent and this trend was also visible for Junb, Plk2, Prune2, and Fbln5 (Figure 5a). To mimic the induction of hyperfiltration, PTECs were pre-exposed to physiological levels of shear (0.25 dyn/cm 2 ) for 4 hr, followed by 16 hr shear stress at the same physiological level or at pathological levels of shear (2.0 dyn/cm 2 ). Expression of Wisp1, Map3k20, Map4k4, Ets1, and Fbln5 was significantly higher at pathological levels of shear compared to physiological levels, while this trend was also visible for Junb and Plk2 (Figure 5b). For the downregulated genes, Stat1 and Jak2, there was no difference in expression between physiological and pathological shear. So, our data indicate that higher levels of shear and a switch from physiological to pathological shear, result in increased gene expression, at least for most the genes analyzed in this experiment.

| DISCUSSION
In this study we used RNA sequencing to get a comprehensive overview of the transcriptome alterations upon fluid shear stress in proximal tubular epithelial cells. Physiological shear stress in renal epithelial cells is ranging from 0.05-1.0 dyn/cm 2 , where proximal tubular cells experience the highest range of shear stress (Essig et al., 2001;Grabias & Konstantopoulos, 2012Weinbaum et al., 2010). We applied a fluid shear stress of 2.0 dyn/cm 2 , which is known to be an increased physio-pathological shear stress, mimicking hyperfiltration after renal mass reduction or during progression of renal disease. Our genome wide RNA sequencing data confirmed previously reported fluid flow-induced changes in gene expression of Cox2 (Ptgs2), Ccl2 (Mcp1), Edn1, Egr1, Snai1, and Cdh1 in renal epithelial cells (Flores et al., 2011(Flores et al., , 2012Maggiorani et al., 2015;Pandit et al., 2015;Schwachtgen et al., 1998). Furthermore, our data reveal >1,500 other genes to be altered by fluid shear stress in PTECs. We validated a subset of genes by qPCR and showed that the shear stress response was time dependent within the first 16 hr. Furthermore, after removal of the shear, the shear-induced gene expression was reversible for some of the genes, while other genes showed similar or higher differential gene expression upon static post incubation. Differences in signaling and cytokine production upon shear may explain the different responses as well as differences in transcriptional activation and stability of transcripts. For example, Fn1 is a very long transcript, which requires more time for transcription and degradation, while Pai1 (Serpine1) has an faster turn-over (Kunnen et al., 2017;'t Hoen et al., 2011).
Pathway analysis indicated increased expression of cell-cell/ cell-matrix interaction genes, including cytoskeletal components, cell adhesion and tight junction molecules, extracellular matrix components and integrins. This suggests strengthening of epithelial cells and their surroundings to resist (increased) physiological shear stress (Duan et al., 2008;Essig et al., 2001;Jang et al., 2013). Another study showed loss of epithelial cell morphology during high pathological shear stress of 5 dyn/cm 2 (Maggiorani et al., 2015). Long-term high shear exposure therefore might also lead to fibrotic deposition and tubulointerstitial lesions, which is commonly seen after renal mass reduction or during FIGURE 2 qPCR validation of RNA sequencing results. Gene expression (log 2 fold change) of selected target genes is altered upon 16 hr fluid shear stress, as measured by quantitative PCR. Parallel plate flow-chamber induced fluid shear stress at 2.0 dyn/cm 2 in PTECs; n = 13 per condition; Hprt served as housekeeping gene to correct for cDNA input; data were normalized to static controls (log 2 fold change = 0). *Indicates significantly altered expression by flow versus no flow (p < 0.05) using a one sample t-test progression of renal diseases (Essig & Friedlander, 2003;Essig et al., 2001;Grabias & Konstantopoulos, 2014;Rohatgi & Flores, 2010;Venkatachalam et al., 2010). Pro-apoptotic as well as pro-survival and cell cycle arrest genes were induced by shear stress, while key players in apoptosis (Bad, Bak, Bax, and caspases) and cell cycle (cyclins and CDKs) were not altered in gene expression. This suggests that apoptosis and cell cycle related gene expression are not dramatically altered during shear exposure.
Core signaling pathways altered by shear stress comprise MAPK and TGF-β signaling. Even more, TGF-β/ALK5-induced target gene expression in renal epithelial cells is partially restrained by MEK1/2mediated signaling (Kunnen et al., 2017). Using ALK4/5/7 inhibitors, we showed that many genes, but not all genes, are dependent on shear induced TGF-β/ALK5 signaling, including genes involved in other core signaling pathways like MAPK and Wnt signaling. The role of TGF-β as a master regulator of the shear response is related to the TGF-β/ALK5 interaction since we previously showed that also TGF-β neutralizing antibodies inhibit the response (Kunnen et al., 2017). It is conceivable that under flow conditions TGF-β processing and binding of the active ligand is enhanced. Interestingly, a recent publication showed that TGF-β can be released from its latency-associated peptide (LAP) by shear stress, probably by forces exerted on α v -β 6 -integrins via the actin cytoskeleton (Dong et al., 2017;Ha, 2017). We also noticed an   , 1999;Lee et al., 2007;Muthusamy et al., 2015), thereby modulating the response to shear. Our data show that the shear stress response of a subset of genes is attenuated upon MEK1/2 inhibition, while other genes showed an enhanced response. Since there are multiple interactions between TGF-β and MAPK/ERK signaling pathways, the integration of these pathways is complex and biological context dependent, and therefore difficult to predict (Kunnen et al., 2017).
In addition to TGF-β signaling, increased expression of other cytokines observed in our study suggests attraction and activation of macrophages and inflammatory cells upon shear in vivo. This is a common phenomenon during development of kidney diseases, where shear stress is fluctuating due to changes in glomerular filtration rate, tubular hyperfiltration and obstruction (Akchurin & Kaskel, 2015).
Altered expression of other growth factors or cytokine signaling pathways include FGF, HB-EGF, PDGF, CXC, and other cytokines.
FGF, HB-EGF, and PDGF can bind to tyrosine kinase receptors that upon activation stimulate the Ras/Raf/ERK (MAPK) pathway and/or the PI3K/ AKT pathway (up-regulated upon shear stress) and/or STAT-signaling (down-regulated upon shear stress) (Pileri & Piccaluga, 2012;Turner & Grose, 2010). At several levels these pathways can be amplified or negatively modulated, and they can interact with each other as well.
Multiple ligand isoforms can bind to the receptors with different affinities. Upon fluid flow, transcript levels of several ligands is increased (Fgf1, Fgf9, Hbegf, Pdgfa, Pdgfb, Pdgfc), but not the receptors. Whether increased signaling is related to endocrine/paracrine loops, as seen for TGF-β, needs more extensive investigation. Interestingly, we also observed altered expression of proteoglycans, like syndecans and glypican, as well as modifying enzymes involved in glycosaminoglycan, heparan-sulphate or chondroitin-sulphate metabolism, which are all involved in glycocalyx remodeling (Reitsma, Slaaf, Vink, van Zandvoort, & oude Egbrink, 2007). Cell-surface-associated heparan sulfate proteoglycans have been shown to be essential for FGF signal transduction and, more general, the glycocalyx is able to significantly modify the cellular response to growth factors including PDGF and FGF. It has been shown that the glycocalyx plays an important role in mechanotransduction of shear stress in endothelial cells. It is required for the cytoskeleton to respond to shear stress and acts as a signaling platform integrating shear stress, growth factor, chemokine and cytokine signaling (Ebong, Lopez-Quintero, Rizzo, Spray, & Tarbell, 2014;Thi, Tarbell, Weinbaum, & Spray, 2004;Zeng, 2017;Zeng & Liu, 2016). So, our data indicate that fluid shear stress induce genes involved in glycocalyx remodeling in PTECs, although it has to be further investigated whether the glycocalyx is equally involved in mechano-sensing upon shear stress in renal epithelial cells.
The shear stress response in PTECs can be regulated by a variety of mechano-sensors at different sub-cellular locations (Curry & Adamson, 2012;Ingber, 2006;Petersen et al., 2016). We investigated the role of primary cilia, since defects in cilia formation and function are associated with developmental disorders and (kidney) diseases (Goetz & Anderson, 2010;Quinlan et al., 2008 Yamamoto & Ando, 2015). This was dependent on AMPK, which is an important kinase in energy metabolism (Carling, 2004;Fisslthaler & Fleming, 2009)  conditions, cells were directly exposed to the indicated level of fluid shear stress. (b) Cells were first pre-exposed for 4 hr to low levels of shear stress (0.25 dyn/cm 2 ), followed by 16 hr shear stress exposure at indicated levels. Expression of all genes was significantly increased by shear stress compared to static controls (dashed line in a) and was flow rate dependent for most genes. (a, b) Parallel plate flow-chamber induced fluid shear stress of 0.25-2.0 dyn/cm 2 in PTECs; t = 16 hr (a) or t = 4 + 16 hr (b); qPCR, Hprt served as housekeeping gene to correct for cDNA input; data normalized to unstimulated controls (fold change); n = 3 per condition. # Significant difference compared to unstimulated control (dashed line in a) or *significant difference between treatment groups (p < 0.05 by one-way ANOVA, followed by post-hoc Fisher's LSD multiple comparison) KUNNEN ET AL.
| 3625 2017). We hypothesize that large variations in shear stress, occurring in kidney diseases, might contribute to the disease phenotype. This hypothesis is supported by our data showing that the expression of several genes involved in TGF-β, MAPK, and Wnt signaling is further elevated upon switching from physiological to pathological levels of shear.
In conclusion, this study provides a