Gene expression in response to cyclic mechanical stretch in primary human dermal fibroblasts.

The human dermal skin is permanently exposed to mechanical stress, for instance during facial expression, which might cause wrinkles with age. Cyclic mechanical stretching of cells results in cellular and cytoskeleton alignment perpendicular to the stretch direction regulating cellular response. With gene expression profiling it was aimed to identify the differentially expressed genes associated with the regulation of the cytoskeleton to investigate the stretch-induced cell alignment mechanism. Here, the transcription activity of the genome in response to cyclic mechanical stress was measured using DNA microarray technology with Agilent SurePrint G3 Human GE 8x60k Microarrays, based on the overall measurement of the mRNA. Gene expression was measured at the beginning of the alignment process showing first reoriented cells after 5 h stretching and at the end after 24 h, where nearly all cells are aligned. Gene expression data of control vs. stretched primary human dermal fibroblasts after 5 h and 24 h demonstrated the regulation of differentially expressed genes associated with metabolism, differentiation and morphology and were deposited at http://www.ncbi.nlm.nih.gov/geo with the accession number GSE58389.

Direct link to deposited data

Deposited data can be found here: http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE58389.

Experimental design, materials and methods

Experimental design

Cell alignment is one of the main cellular responses to cyclic mechanical uniaxial stretch [1] and might be associated with mechanically induced wrinkle formation in the skin [2]. The identification of the mechanism leading to alignment may allow analysis and modulation of its role in the formation of mechanically induced wrinkles. To identify changes in gene expression associated with mechanical stretch-induced cell alignment, a whole genome microarray study was performed on primary human dermal fibroblasts (PHDF) subjected to cyclic uniaxial stretching using a Flexer®Cell Tension Plus System. Gene expression was measured at the beginning of the alignment process showing first reoriented cells after 5 h stretching and at the end after 24 h, when nearly all cells are aligned perpendicular to the stretch-direction (Fig. 1). In total, PHDF from ten donors were cultured on BioFlex culture plates and stretched for 5 h and 24 h or left untreated as controls to account for changes occurring during cell culture. This resulted in 4 samples for each of the subjects (control/treated and 5 h/24 h), i.e. 40 samples in total (Table 1).

Fig. 1

PHDF were cultured on BioFlex culture plates for 24 h and then stretched with the FX-4000T™ Tension Plus™ System. Cell morphology of PHDF in response to 24 h cyclic mechanical stretching was observed with transmitted light microscopy. Micrographs of static cells show randomized orientated PHDF. In contrary, stretched cells were oriented perpendicular to stretch direction.

Table 1

Experimental design of microarray-based gene expression profiling. PHDF from ten subjects were cultured on BioFlex culture plates and stretched for 5 h and 24 h or left untreated. RNA was isolated and used for gene expression profiling. Each subject provided the 4 samples with the control/treatment combinations (control; stretched).

	Subject ID	Primary human dermal fibroblasts isolated from 10 subjects
Control		Stretched
Sample name	Array ID	Sample name	Array ID
5 h	392	Reu_1	co_5h_392	Reu_21	tr_5h_392
5 h	420	Reu_3	co_5h_420	Reu_23	tr_5h_420
5 h	464	Reu_5	co_5h_464	Reu_25	tr_5h_464
5 h	425	Reu_7	co_5h_425	Reu_27	tr_5h_425
5 h	432	Reu_9	co_5h_432	Reu_29	tr_5h_432
5 h	416	Reu_11	co_5h_416	Reu_31	tr_5h_416
5 h	465	Reu_13	co_5h_465	Reu_33	tr_5h_465
5 h	446	Reu_15	co_5h_446	Reu_35	tr_5h_446
5 h	387	Reu_17	co_5h_387	Reu_37	tr_5h_387
5 h	445	Reu_19	co_5h_445	Reu_39	tr_5h_445
24 h	392	Reu_2	co_24h_392	Reu_22	tr_24h_392
24 h	420	Reu_4	co_24h_420	Reu_24	tr_24h_420
24 h	464	Reu_6	co_24h_464	Reu_26	tr_24h_464
24 h	425	Reu_8	co_24h_425	Reu_28	tr_24h_425
24 h	432	Reu_10	co_24h_432	Reu_30	tr_24h_432
24 h	416	Reu_12	co_24h_416	Reu_32	tr_24h_416
24 h	465	Reu_14	co_24h_465	Reu_34	tr_24h_465
24 h	446	Reu_16	co_24h_446	Reu_36	tr_24h_446
24 h	387	Reu_18	co_24h_387	Reu_38	tr_24h_387
24 h	445	Reu_20	co_24h_445	Reu_40	tr_24h_445

Materials and methods

Cell culture

Dermal fibroblasts were isolated from skin biopsies obtained from plastic surgery. The biopsies were cut in stripes and incubated in dispase II for 3 h at 37 °C. Afterwards the dermis was separated from the epidermis, cut into smaller pieces and digested in collagenase I o/n at 37 °C. The suspension was filtered and the primary human dermal fibroblasts in the filtrate were seeded in cell culture flasks. For long term storage the PHDF were cultured until passage one in cell culture flasks, harvested by centrifugation (5 min, 1000 g) after Trypsin/EDTA digestion for 5 min at 37 °C and resolved in freezing medium. PHDF were thawed with prewarmed DMEM containing 10% FCS, 2 mM GlutaMax™-I and 0.1 mg/ml penicillin/streptomycin, plated in cell culture flasks and incubated at 37 °C in an atmosphere of 5% CO2 and 90% humidity until confluence. Then cells were subcultured in BioFlex culture plates.

Cyclic mechanical stretching in vitro

Mechanical stretching was performed on flexible silicon membranes using FX-4000T™ Tension Plus™ System. PHDF were seeded on BioFlex culture plates coated with Collagen I and cultured until subconfluence of 70% and cells were serum deprived o/n. For stretching experiments culture plates were mounted on the Baseplate™. Cyclic stretch was applied in the FX-4000T™ Tension Plus™ System with 16% elongation, 0.5 Hz in a half sinus regimen. By application of a vacuum a depression occurs and the silicon membranes with adhering cells were stretched over the loading posts. Cell alignment was microscopically observed at the outer circular region of the well. At 5 h and 24 h cells were harvested for mRNA isolation.

RNA isolation from primary human dermal fibroblasts

PHDF from 10 donors were cultured on BioFlex culture plates for 24 h and stretched for 5 h and 24 h or left untreated. To separate the inhomogeneous stretching areas of the BioFlex culture plates the silicon membranes were punched with a 2 cm diameter punch. Isolation of RNA was done with RNeasy Mini Kit according to manufacturer product information from the outer circular region of the well. DNA and RNA were precipitated with 70% ethanol and bound to a silica membrane. DNA was digested using DNAse I. RNA was eluted with 30 μl RNAse free water and subjected to Experion automated electrophoresis for quality control.

Microarray hybridization and data processing

The obtained RNA was transcribed into cDNA and than subjected to microarray hybridization. 100 ng of each total RNA sample was used for the linear T7-based amplification step. To produce Cy3-labeled cRNA, the RNA samples were amplified and labeled using the Agilent Low Input Quick Amp Labeling Kit (Agilent Technologies) following the manufacturer's protocol. Yields of cRNA and the dye-incorporation rate were measured with ND-1000 Spectrophotometer (NanoDrop Technologies). The hybridization procedure was performed according to the Agilent 60-mer oligo microarray processing protocol using the Agilent Gene Expression Hybridization Kit (Agilent Technologies) (Design ID 028004). Subsequently, 600 ng Cy3-labeled fragmented cRNA in hybridization buffer was hybridized 17 h at 65 °C to Agilent SurePrint G3 Human GE 8x60k Microarrays using Agilent's recommended hybridization chamber and oven. Finally, the microarrays were washed once with the Agilent Gene Expression wash buffer 1 for 1 min at room temperature followed by a second wash with preheated Agilent Gene Expression wash buffer 2 (37 °C) for 1 min. The last washing step was performed with acetonitrile. Fluorescence signals of the hybridized Agilent microarrays were detected using Agilent's Microarray Scanner System (G2505C, Agilent Technologies, Palo Alto, USA). The Agilent Feature Extraction Software (FES 10.7.3.1) was used to obtain and process the microarray image files. The preprocessing started with the conversion of the data after using Agilent Feature Extraction software in txt files suitable for all subsequent analysis steps, which were mainly performed with the R statistical software and its Bioconductor packages. For the preprocessing the agi4x44kpreprocess package was used [3]. To this end an annotation package for the Agilent Whole Genome 8x60k chip has previously been created. The annotation package for the 8x60k Agilent chip was created using Bioconductor SQLForge and AnnotationDbi packages based on annotation information provided as an Excel file by Agilent. After loading the data files into R, all data were background-corrected. Data were then normalized between arrays using quantile-normalization and transformed to log2-scale, which enabled comparison of samples loaded on different arrays. After normalization, a set of quality control steps was performed to filter low-quality probes. Filtering of probes was based on quality flags set by the Agilent Feature Extraction Software. Most of the probes (~ 42%) were filtered because of not being sufficiently above background.

Between sample comparisons of whole genome expression profiles

The processed and filtered data were subjected to between sample comparisons. To evaluate consistency for batch processing, between sample Euclidian distances were calculated using all expression values, scaled by row, and plotted as a false-color heat map with red indicating high differences and blue indicating low differences between compared samples (Fig. 2). The heat map showed lower differences between samples from the same subject, and otherwise relatively homogenous coloring demonstrating that gene expression of samples differed to similar extents indicating no batch effects. Differentially expressed genes were identified using the empirical Bayes method implemented in the limma package for the R statistical software [4], [5]. Paired statistical analysis of control and treated cells at 5 h and 24 h was performed with a Benjamini–Hochberg adjusted p-value cut-off of ≤ 0.05 [6], [7], [8]. Significantly expressed gene probes were filtered to exclude gene probes with a log(fold change) ≤ 0.5 to improve separation resulting in 603 differentially expressed gene probes for 5 h treatment and 177 differentially expressed gene probes for 24 h treatment (for 24 h Table 2, for 5 h Supplementary data, not assignable gene probes were excluded). Heat maps were drawn for significantly differentially expressed gene probes for 5 h treatment versus control (Supplementary data) and 24 h treatment versus control comparisons (Fig. 3), using gene expression values scaled by row, with red color indicating upregulation relative to control and blue color indicating downregulation with respect to control. The heat map of significantly differentially expressed gene probes for 5 h shows a less clear separation of samples according to their treatment, control or stretched, indicating that gene expression changes in response to 5 h stretching were less prominent. It can be assumed that this is in correlation with the observed morphological phenotype (Fig. 1). The cell alignment process starts at 5 h and therefore the cells are still in the reorientation process at this time point and gene expression changes are unassertive. Heat map of differentially expressed gene probes for 24 h demonstrates a clear separation of samples according to their treatment indicating that after 24 h stretching gene expression was explicitly changed at that time point, when cell alignment process was completed.

Fig. 2

Between sample comparisons of gene expression changes measured with Whole genome 8x60k array chip hybridization. PHDF from 10 donors were cultured on BioFlex culture plates and then stretched for 5 h and 24 h or left untreated. RNA was isolated and subjected to microarray-based whole genome expression profiling. Fluorescence intensities were preprocessed with statistical software R using its Bioconductor package. The false-color heat map was drawn in R representing clustering of samples by the means of between sample distances. The scaled expression value, denoted as the row Z-score, is plotted in dark red–dark blue color code.

Table 2

Table of significant differentially expressed genes at 24 h with scaled log(fold changes) and p-values. Paired statistical analysis of control and treated cells at 24 h was performed with a Benjamini–Hochberg adjusted p-value cut-off of ≤ 0.05. Significantly expressed gene probes were filtered to exclude gene probes with a log(fold change) ≤ 0.5 to improve separation.

24 h treatment (115 gene probes)
Log(fold changes)	p-Value	Description
0.5398	0.000671083	MIR143 host gene (non-protein coding)
0.5328	1.29E− 05	Uncharacterized LOC100506860
− 0 .7350	6.04E− 06	Uncharacterized LOC100506377
0 .5595	1.11E− 06	Long intergenic non-protein coding RNA 467
0 .5431	3.37E− 09	Uncharacterized LOC440993
0 .5298	1.15E− 07	Uncharacterized LOC100499194
0 .6250	0.000297072	Uncharacterized LOC100505687
0.9110	2.44E− 05	Kin of IRRE like 3 (Drosophila)
− 0 .7072	3.16E− 07	Niemann –Pick disease, type C1
− 0.7884	2.04E− 07	Methylsterol monooxygenase 1
0 .5714	6.93E− 06	Cysteine-rich protein 2
0 .5357	2.96E− 06	Malic enzyme 3, NADP(+)-dependent, mitochondrial
0 .5403	1.90E− 05	Regulator of G-protein signaling 11
− 0 .6491	5.66E− 09	Hydroxysteroid (17-beta) dehydrogenase 7
0 .5137	0.000527329	Kruppel-like factor 2 (lung)
− 0.5687	1.85E− 05	Tumor necrosis factor receptor superfamily, member 14
− 0 .6834	3.44E− 09	Syntaxin 3
− 0 .5257	1.85E− 08	Family with sequence similarity 134, member A
0 .5115	0.000514841	Slit homolog 2 (Drosophila)
− 0.7784	3.27E− 08	Glucosamine-6-phosphate deaminase 1
− 0 .5701	4.65E− 07	Squalene epoxidase
− 0 .9861	0.000815157	Ferritin, heavy polypeptide-like 17
0.6104	6.10E− 05	Synaptogyrin 4
− 0 .5113	3.03E− 07	5 ′–3′ exoribonuclease 2
− 0.5325	9.73E− 10	Mannose-6-phosphate receptor (cation dependent)
0 .7487	1.92E− 11	AHNAK nucleoprotein
0 .7842	7.31E− 05	Early growth response 1
0 .5585	0.000882563	Tumor necrosis factor receptor superfamily, member 6b, decoy
0.6287	9.97E− 11	Plectin
0 .5429	1.26E− 07	MAD2 mitotic arrest deficient-like 2 (yeast)
0 .6556	2.76E− 05	Chromosome 10 open reading frame 54
0 .5536	0.000577494	Cell division cycle associated 7
0.5976	2.14E− 05	Creatine kinase, brain
− 0 .5909	3.60E− 05	Chromosome 14 open reading frame 1
0 .5561	0.000729642	Rhomboid, veinlet-like 1 (Drosophila)
0.5233	0.000842352	Late cornified envelope 1C
− 0.5510	6.57E− 06	3-Hydroxy-3-methylglutaryl-CoA reductase
− 0 .5259	1.02E-07	Basic helix–loop–helix domain containing, class B, 9
− 0.6332	8.54E-07	ATPase, H + transporting, lysosomal 56/58 kDa, V1 subunit B2
− 0.5538	4.46E-05	StAR-related lipid transfer (START) domain containing 4
− 0.6156	1.36E-07	Transmembrane protein 217
− 0.5978	4.37E-05	Uncharacterized LOC197187
0.5690	0.000120951	Chromosome 14 open reading frame 80
0.5223	4.57E-05	Interleukin 17 receptor E
0.5694	0.000210521	Tubulin, alpha 3d
− 0.5380	4.86E-05	Sec23 homolog B (S. cerevisiae)
0.5433	1.87E-05	MYCN opposite strand/antisense RNA (non-protein coding)
− 0.5530	4.25E-06	Malic enzyme 1, NADP(+)-dependent, cytosolic
0.5807	1.85E-05	Thioredoxin interacting protein
− 0.5730	1.20E-05	Multiple coagulation factor deficiency 2
− 0.6702	7.29E-08	Isopentenyl-diphosphate delta isomerase 1
0.6304	0.000232544	Ubiquitin-like with PHD and ring finger domains 1
− 0.5725	1.55E-08	Ankyrin repeat domain 37
− 0.5130	3.22E-07	Sorting nexin 3
0.6085	0.000458856	Fibrosin-like 1
− 0.6386	6.17E-05	RAB9B, member RAS oncogene family pseudogene 1
− 0.5078	1.77E-05	dCMP deaminase
0.7674	0.000490824	Ribonucleoprotein, PTB-binding 1
− 0.6215	4.02E-08	Sialidase 1 (lysosomal sialidase)
− 0.6531	7.17E-06	Solute carrier family 43, member 3
0.5878	0.000128578	Chromosome 19 open reading frame 55
− 0.5231	0.000172375	3-Hydroxy-3-methylglutaryl-CoA synthase 1 (soluble)
0.5502	0.000793717	Chromosome 1 open reading frame 229
− 0.5190	0.000100862	Junctional adhesion molecule 3
0.6320	0.000533996	Histone cluster 1, H3i
− 0.6842	0.000136983	Proprotein convertase subtilisin/kexin type 9
− 0.5162	8.24E− 05	NECAP endocytosis associated 1
− 0 .7197	3.99E− 06	Stearoyl-CoA desaturase (delta-9-desaturase)
− 0 .5764	2.22E− 07	Wilms tumor 1 associated protein
0 .5780	0.000785569	Histone cluster 1, H2am
0.5660	1.36E− 06	Long intergenic non-protein coding RNA 577
0 .5647	0.000804159	Secretoglobin, family 3A, member 1
0.5104	1.32E− 05	PDZ and LIM domain 7 (enigma)
− 0 .5876	8.43E− 05	Insulin induced gene 1
0 .5709	0.000973297	LIM homeobox 3
0.5330	0.000551954	v-Raf murine sarcoma 3611 viral oncogene homolog 2, pseudogene
− 0.5549	2.69E− 07	Hydroxysteroid (17-beta) dehydrogenase 7
− 0 .5286	1.51E− 09	NAD(P) dependent steroid dehydrogenase-like
− 0 .5962	7.70E− 05	Paraneoplastic Ma antigen 1
− 0 .6144	0.000269973	Potassium voltage-gated channel, Isk-related family, member 4
− 0.5522	6.60E− 10	UBX domain protein 8
− 0 .6874	0.000809869	Small nucleolar RNA, C/D box 3B-1
− 0.5275	1.04E− 07	Patched domain containing 3 pseudogene
0 .7635	0.001070526	Histone cluster 2, H3a
0.5250	0.001093476	NCK adaptor protein 2
0.5829	0.000126935	Uncharacterized protein FLJ25694
0.5291	0.001058371	Homeobox A10
0.5255	2.84E− 05	Collagen, type VIII, alpha 2
0 .5281	0.000314357	Uncharacterized LOC643549
0.5306	0.000422623	Calcyon neuron-specific vesicular protein
0.5472	1.86E− 05	Hairy and enhancer of split 6 (Drosophila)
− 0 .5958	0.000268984	Fas (TNF receptor superfamily, member 6)
− 0.5157	1.91E− 07	Heat shock 105 kDa/110 kDa protein 1
0 .5736	0.00096652	Keratin 7
0.5345	1.61E− 10	Chromosome 1 open reading frame 86
− 0 .6078	1.97E− 07	Activin A receptor, type I
0 .5641	0.000259433	High mobility group protein HMG-I/HMG-Y-like
− 0.5045	7.21E− 12	Putative homeodomain transcription factor 1
0 .5329	5.05E− 05	Uncharacterized LOC388152
0 .5787	3.07E− 05	ADAMTS-like 4
0 .6463	6.78E− 05	Exocyst complex component 7
0 .5300	2.30E− 07	Exocyst complex component 7
0 .5380	0.000583772	Transmembrane protein 8C
0.5741	5.36E− 05	Protein kinase, membrane associated tyrosine/threonine 1
0 .5195	0.000353071	Golgin A6 family-like 1 pseudogene
− 1.0093	0.000147492	Ubiquitin specific peptidase 6 (Tre-2 oncogene)
0.5359	0.000519288	NLR family, CARD domain containing 3
0.6446	0.000549649	CD72 molecule
0.5893	2.66E− 05	B1 for mucin
0 .6576	1.74E− 05	Small nucleolar RNA, H/ACA box 71A
− 0 .5154	1.08E− 08	Cysteine and histidine-rich domain (CHORD) containing 1
0 .5014	1.23E− 05	Uncharacterized LOC100507637
− 0 .6767	3.14E− 06	Chloride intracellular channel 2
0.5141	0.001075593	Polycystic kidney disease 1 (autosomal dominant)
0.5094	0.000108073	Protection of telomeres 1 homolog (S. pombe)

Fig. 3

The scaled gene expression for control versus stretched samples at 24 h of each probe is shown as the row Z-score; it is plotted in a red–blue color scale with red indicating high expression and blue indicating low expression. Significantly differentially expressed gene probes were filtered to exclude gene probes with a log(fold change) ≤ 0.5. Heat map of differentially expressed gene probes for 24 h  demonstrates a clear separation of samples according to their treatment.

Specifications
Organism/cell line/tissue	Primary human dermal fibroblasts isolated from skin biopsies
Sequencer or array type	Agilent SurePrint G3 Human GE 8x60k Microarray (G4858A-028004)
Data format	Normalized signal intensity
Experimental factors	Static vs. stretched cells
Experimental features	Comparison of gene expression profiles of cyclic mechanical stretched primary human dermal fibroblasts with untreated control cells isolated from ten donors after 5 h and 24 h stretching
Consent	Donors provided written consent