Next-generation sequencing of small RNAs from inner ear sensory epithelium identifies microRNAs and defines regulatory pathways

Background The mammalian inner ear contains sensory organs, the organ of Corti in the cochlea and cristae and maculae in the vestibule, with each comprised of patterned sensory epithelia that are responsible for hearing and balance. The development, cell fate, patterning, and innervation of both the sensory and nonsensory regions of the inner ear are governed by tight regulation involving, among others, transcription factors and microRNAs (miRNAs). In humans, mutations in specific miRNA genes are associated with hearing loss. In mice, experimental reduction or mutations of miRNAs in the inner ear leads to severe developmental and structural abnormalities. A comprehensive identification of miRNAs in the sensory epithelia and their gene targets will enable pathways of auditory and vestibular function to be defined. Results In this study, we used Next-Generation Sequencing (NGS) to identify the most prominent miRNAs in the inner ear and to define miRNA-target pairs that form pathways crucial for the function of the sensory epithelial cells. NGS of RNA from inner ear sensory epithelial cells led to the identification of 455 miRNAs in both cochlear and vestibular sensory epithelium, with 30 and 44 miRNAs found in only cochlea or vestibule, respectively. miR-6715-3p and miR-6715-5p were defined for the first time in the inner ear. Gene targets were identified for each of these miRNAs, including Arhgap12, a GTPase activating protein, for miR-6715-3p, implicating this miRNA in sensory hair cell bundle development, actin reorganization, cell adhesion and inner ear morphogenesis. Conclusions This study provides a comprehensive atlas of miRNAs in the inner ear sensory epithelia. The results provide further support of the essential regulatory role of miRNAs in inner ear sensory epithelia and in regulating pathways that define development and growth of these cells. Electronic supplementary material The online version of this article (doi:10.1186/1471-2164-15-484) contains supplementary material, which is available to authorized users.

Background miRNAs play an essential role in inner ear development [1]. miRNAs are small non-coding RNAs that regulate gene expression post-transcriptionally. By binding to sequences in the 3' untranslated region (UTR) of mRNAs, a miRNA can inhibit target mRNAs by translational suppression and mRNA destabilization [2]. miRNAs have been implicated in hearing loss in humans, since mutations in the seed region of miR-96 or its gene are associated with hearing loss in three extended families [3,4]. miRNAs were also reported in other human ear pathologies, including an elevation of miR-21 in cholesteatomas [5] and in vestibular schwannomas [6]. miRNAs were found to regulate the otitis media inflammatory response [7]. Their study in humans, however, has been hampered by the unavailability of inner ear RNA from human subjects, making the mouse an invaluable model for studying miRNA development and regulation in the inner ear [8]. In an ENUinduced mutant mouse exhibiting deafness and vestibular dysfunction, a mutation in the miR-96 seed region was found to be responsible for hearing loss, presumably due to the alteration in expression and hence function of both direct and indirect gene targets [9]. A number of mouse mutants affecting miRNA regulation through Dicer have been created, leading to developmental inner ear defects [10,11]. These conditional knock-out mouse mutants have been instrumental in demonstrating that miRNAs are vital for inner ear morphogenesis and development of the sensory epithelia and sensory neurons.
Expression of over 300 miRNAs have been reported in mouse and rat inner ears thus far, identified by microarray analysis [8,10,12,13], with the majority of work performed on the triad mir-96, −182, and −183 [14]. In order to further identify and characterize miRNAs in the mammalian inner ear, we used NGS for the first time to identify miRNAs in cochlear and vestibular sensory epithelia. Using this method, we identified over 500 inner ear miRNAs, both in common between and unique to each tissue. We further validated and characterized the expression of two miRNAs, miR-6715-3p and miR-6715-5p, in the mouse inner ear. Vezatin, an adherens junctions transmembrane protein previously implicated in deafness [15] and Arhgap12, a Rho GTPase activating protein [16], were verified to be gene targets of miR-6715-3p. Protocadherin 19, a member of the δ2 subclass of nonclustered protocadherins and associated with epilepsy and mental retardation [17], is a gene target of miR-6715-5p. Arhgap12, newly described in the inner ear, is expressed in the epithelial cell-cell junctions of the hair and supporting cells. The miR-6715-3p-Arhgap12 miRNA target pair defines a new regulatory pathway. Understanding these interactions may shed light on known and novel miRNAs, their effect on the development of normal and impaired hearing, and the mechanisms leading towards deafness.

RNA-seq derived miRNA transcription profile of the mouse inner ear sensory epithelium
In order to identify miRNAs in the auditory and vestibular sensory epithelium, we sequenced short RNA molecules from cochlear and vestibular sensory epithelia by high-throughput RNA sequencing (RNA-seq). For the cochlea, the sensory epithelium consisted of the organ of Corti, including hair cells, supporting cells, and cells of the greater and lesser epithelial ridges, as well as mesenchymal and neuronal cells. For the vestibule, the sensory epithelium was derived from the saccule, utricle and the lateral and anterior cristae and included hair cells and supporting cells. The advantage of using high-throughput sequencing is that this technique allows for assessment of transcript expression in a given sample at a large dynamic range, from minute quantities to highly abundant ones, and does not require a-priori knowledge for the identification of modified and novel RNA molecules in addition to known ones. cDNA libraries from small RNAs were prepared from dissected sensory epithelium of post-natal day (P)0 mice. Subsequent analysis of sequence reads was done using miRNAkey, a tool for miRNA differential expression analysis [18]. Processed reads were aligned to the mature Mus musculus miRNA database (http://www.mir base.org), allowing one mismatch to occur between a read and the reference, due to isomiRs [19]. Sixty-six and 61.4% of the total reads from cochlear and vestibular sensory epithelia, respectively, were aligned to known mouse miRNAs, indicating a successful isolation of miRNAs. The majority of mapped reads in cochlear and vestibular sensory epithelia, 96.9% in the cochlea and 97.6% in the vestibule, were uniquely mapped to a single miRNA. Overall, 485 and 499 known miRNAs had an RPM (reads per million) mapping higher than 1 in the cochlear and vestibular sensory epithelium, respectively (Additional file 1: Table S1). These lists include all miRNAs that have been studied thus far in the mouse inner ear sensory epithelia, as well as miRNAs that have not been previously detected. Nine miRNAs were most highly expressed in the cochlea and the vestibule, with a similar order in both, except for miR-204-5p, which was the fourth most prevalent in the cochlea and the lowest in the vestibule (Table 1).
Differential expression among the known miRNAs was inferred by miRNAkey, based on the widely-accepted assumption that the number of reads mapped to each miRNA is relative to its level of expression [20,21]. The reads numbers in each sample were normalized to RPM values, and the relative expression levels of the different miRNAs in the two samples were calculated. In general, miRNA expression was similar in the cochlear and vestibular sensory epithelia. The majority of miRNAs with RPM > 1 were found in both tissues (455, or 86%), with 30 found only in the cochlea and 44 found only in the vestibule ( Figure 1A, Additional file 1: Table S1). The rates of the highest expressing miRNAs in each tissue were also highly similar between tissues (Pearson coefficient = 0.99; P < 10 −8 ) ( Figure 1B). We defined a miRNA as differentially expressed if the fold change between samples was higher than 4 and P < 10 −4 (Bonferroni corrected). Eleven miRNAs were differentially expressed in cochlear and vestibular sensory epithelia ( Table 2).
Targets for most highly expressed miRNAs miRNAs and targets exert their function through regulation of biological processes. As each miRNA can have multiple targets, we performed an analysis to discover additional mRNA genes that are computationally predicted to be targets of the three most highly expressed miRNAs, miR-182, miR-181a and miR-26a. Among numerous predictions, a few were of particular interest due to their involvement in auditory function. The targets were defined according to their presence in the list of genes that are associated with deafness in the Mouse Genome Informatics (MGI) (http://www.informatics. jax.org/diseasePortal/) database and predicted by two prediction algorithms, TargetScanMouse and miRanda ( Table 3). All miRNA targets presented are involved in auditory function and hearing loss. As a few of the predicted targets overlap, we postulated that there might be underlying regulatory networks prompted by these miR-NAs. This was further illustrated using gene ontology (GO) of biological-based processes, performed by gene-MANIA (http://www.genemania.org), where 6/16 genes fell in the categories 'sensory perception of sound' and 'sensory perception of mechanical stimulus, ' while 5/16 genes fell in the categories 'inner ear development' and 'ear development' (Additional file 2: Table S2). These predictions may set the basis for signaling pathways in the inner ear, subject to experimental validation.

Identification of inner ear miRNAs
In order to identify inner ear miRNAs, the genes associated with "impaired hearing" (MP:0006325) and "deafness" (MP:0001967) were downloaded from MGI (Additional file 3: Table S3). Since these genes are crucial for inner ear development and function, we predicted that miRNAs that reside in their introns may have a role in the inner ear. Out of the hearing associated genes, ten had pri-miRNAs  in their introns (Additional file 4: Table S4). One pri-miRNA, expressed in both cochlear and vestibular sensory epithelia with an RPM > 1, was found in the intron of deafness-related gene β-tectorin (Tectb) (Figure 2A), which is expressed in the mouse inner ear and associated with altered hearing [22]. This pri-miRNA produces two mature miRNAs, miR-6715-3p and miR-6715-5p ( Figure 2B). These miRs were supported by 2.08 and 5.06 reads per million library reads, respectively ( Figure 2C). These miRs were also found in an analysis of sequences from mouse data sets derived from different tissues [23]. Provisional ID chr19_45771, available on miRBase (miR-6715; MI0025026), was detected by an extremely low number of reads (0.0021-0.0126 per million library reads) in oocytes, embryonic day (E)7.5 embryo, E15.5 brain, hippocampus, heart, skin, muscle cells, lung and liver by RNA-seq.

Expression of new inner ear miRNAs
The expression of miR-6715a-3p and miR-6715a-5p, detected by RNA-seq, was validated by qRT-PCR analysis ( Figure 3A). We found both miRNAs to be expressed in the sensory epithelium of the mouse cochlea, at E16, P0 and P8. The expression of both miRNAs was dynamic and increased with age, showing higher expression at P8 compared to P0 and E16, and at P0 compared to E16 (P < 0.005). The two miRNAs were undetected in brain, liver, kidney and heart. This is consistent with the previous RNA-seq detection [23]. The inner ear expression, and the absence or very low expression of the miRNAs in other tested tissues, suggests that they are inner ear-specific miRNAs. The expression prior to the onset of hearing suggests a developmental role for these miRNAs.
In order to explore the spatial expression pattern of the two inner ear miRNAs, in situ hybridization analysis was performed on inner ears derived from P0 C57/BL6 mice, at the same age the RNA-seq was performed. Consistent with the qRT-PCR results, expression of the two miRNAs was observed in P0 mice ( Figure 3B). A clear staining of both miRNAs was observed in many parts of the inner ear, in spiral and vestibular ganglia, basilar membrane and Reissner's membrane, the stria vascularis and spiral ligament. In a higher magnification, in the cochlea, expression of both miRNAs was observed in the organ of Corti, in the spiral limbus and in particular in the hair and supporting cells. In the vestibular system, staining was observed in sensory epithelia of the saccule and the utricle, in hair and supporting cells.

Target pathways in auditory and vestibular function
In order to determine what pathway might be impacted by miR-6715-3p and miR-6715-5p regulation, we searched for potential targets. Targets were predicted by TargetScan Custom, which searches for a complementary 3'UTR given a putative seed sequence. miR-6715-3p and miR-6715-5p were each predicted to target 201 and 206 targets, respectively, with ten overlapping targets (Additional file 5: Table S5). We examined which potential targets were inner ear-expressing genes, genes associated with hearing loss in human or mice, or genes encoding proteins involved in pathways relevant for inner ear function. These included protocadherin 19 (Pcdh19) [24] as a candidate target for  miR-6715-5p, and vezatin (Vezt) [15], Scn8a [25] and Arhgap12 [16] as candidate targets for miR-6715-3p.
To test whether these potential targets have a direct interaction with miR-6715-3p and miR-6715-5p, we examined luciferase activity by inserting the 3'UTR of the potential targets downstream of a luciferase gene. Transfections of these constructs, together with a pre miR-6715 expressing vector into HEK293T cells, and measurement of the relative luciferase activity compared to a mutant 3'UTR construct, were indicative of targeting. Pcdh19 was found to be a direct target of miR-6715a-5p ( Figure 4A), with a reduction of luciferase activity by approximately 35%. Vezt and Arhgap12 were found to be direct targets of miR-6715-3p, with a reduction of luciferase activity by 55% and 15%, respectively ( Figure 4B, C). Snc8a is not a direct target of miR-6715a-3p, as luciferase activity was not reduced ( Figure 4D).
Arhgap12 was chosen for further analysis, due to its putative role as a GTPase activating protein (GAPs) for small G proteins and expression in epithelial cell-cell junctions [16]. Temporal expression of Arhgap12 mRNA was evaluated by qRT-PCR in the cochlear sensory epithelium. In the cochlea, Arhgap12 was expressed at E16, P0 and P8, with a small, though not significant, increase in age ( Figure 5). This protein was not previously characterized in the inner ear. We found Arhgap12 expression, using a commercial antibody examined for specificity by peptide competition (Additional file 6) in the organ of Corti, in hair and supporting cells, the stria vascularis, spiral ligament and Reissner's membrane ( Figure 6A, B). The strongest expression was observed in nerve cells, in the spiral ganglia ( Figure 6C). In the vestibule, expression was found in the hair and supporting cells ( Figure 6D). In whole mount preparations, Arh-gap12 was found to be present in hair and supporting cells ( Figure 7A-C). Its expression was observed in the same hair cell and cuticular plate planes as ZO-1, but did not overlap ( Figure 7A-C). The outer hair cell and Figure 2 The genomic location, predicted structure and relative level of expression of miR-6715-3p and miR-6715-5p. A. pre-miR6715 is located in the second intron of the Tectb gene. B. miRDeep2 prediction of the RNA secondary hairpin structure of the unprocessed miR-6715-3p, miR-6715-5p, and loop. The seed region of miR-6715-3p is indicated in blue and of miR-6715-5p in green. C. Middle density plot shows the distribution of reads in the predicted precursor sequence, as produced by miRDeep2.

Discussion
The function of miRNAs and their targets in the mammalian inner ear are now being discovered, as the link between miRNAs and hearing and deafness has been established. The inner ear is a complex tissue with multiple cell types, with miRNAs involved in fine-tuning the multitude pathways required for functional hearing. While there are a relatively small number of miRNAs in mammals, their role in regulation far surpasses their number, since each miRNA can exert its effect on hundreds of downstream targets [26]. Hundreds of miRNAs have been discovered in the inner ear; only a subset of these miRNAs have validated targets (reviewed in Table 1 [27]).
Our study was designed to profile the general expression of miRNAs in the inner ear and a comparison of the auditory and vestibular epithelia. We chose the time point of p0 since work by our group [10,13] and others have demonstrated that miRNAs are expressed and functional at this stage, including miRNA-96, which is essential for hearing [3,9]. RNA-seq analysis led to the identification of over 500 miRNAs, with miR-182 as the most highly expressed miRNA in both sensory epithelia (Table 1) and accounting for more than 50% of the 20 most highly expressed miRNAs. A few targets have been confirmed for miR-182, including Sox2, Clic5 and Tbx1 [28][29][30]. The second most abundant miRNA in the cochlear and vestibular sensory epithelia, miR-181a-5p, has no known validated targets in the inner ear. This miRNA has been studied in other systems and found to have a proproliferative role in cultured human myeloid leukemia cells [31] and regulates thymic selection in mouse T-cells Figure 3 miR-6715-3p and miR-6715-5p expression in the inner ear. A. qRT-PCR analysis and validation of novel miRNAs in the mouse cochlea sensory epithelium. miR-6715-3p and miR-6715-5p exhibited increased expression with age between E16 and P8 (*) P < 0.05, (**) P < 0.005. n = 5. B. In situ hybridization of P0 mouse cochlear and vestibular sections. miR-6715-3p and miR-6715-5p expression was found in the hair and supporting cells of both cochlea and vestibule, spiral and vestibular ganglia, stria vascularis, basilar membrane and in Reissner's membrane. sa, saccule; ut, utricle; sg, spiral ganglia; sm, scala media; st, scala tympani; sv, scala vestibule; ohc, outer hair cells; ihc, inner hair cells; hc, hair cells. Bar: 100 μm. [32]. The role of this miRNA has been studied in the avian inner ear, where it has been shown to promote proliferation [33].
We validated a number of targets for miR-6715-3p and miR-6715-5p. Predicted targets for the miRNAs described in this study have direct or potential implications for hearing loss. Pcdh19 is expressed in the zebrafish inner ear [24]. Mutations in protocadherin 19 (PCDH19) have been associated with X chromosome-linked epilepsy and mental retardation and it is a member of the protocadherin family of proteins [17]. This family includes protocadherin 15 (Pcdh15), with a known function in the tip links of the stereocilia of the inner ear in conjunction with cadherin 23 [34] and has been associated with multiple mutations in both non-syndromic and syndromic hearing loss in humans [35,36].
Vezatin (Vezt) is an adherens junctions transmembrane protein, and while it is ubiquitously expressed in these junctions, it is especially prevalent in in hair and supporting cells of the mouse inner ear [15]. Vezt mutant mice  suffer from late onset progressive hearing loss, due to hair cell apoptosis. Most compelling, the resistance of the organ of Corti to mechanical stress was reduced when vezatin was absent in the inner ears of conditional mutants. These studies led to the prediction that vezation is required to maintain junction integrity and is crucial to protect the junctions from mechanical stress due to sound trauma.
Arhgap12, a Rho GTPase activating protein (GAP) family member, is a component of cell-cell junctions in kidney epithelial cells [16]. While not previously known to be expressed in the inner ear, its target, Rac1, localizes to hair and supporting cells, specifically in the stereocilia and kinocilium of the hair cells. It was found to be involved in hair cell morphogenesis, specifically in hair bundle development, epithelia actin assembly and cell adhesion. Inactivation of Rac1 specifically in the otic epithelium leads to defects in cochlear morphogenesis, hair bundle orientation and cellular organization of the organ of Corti [37]. Arhgap12 was found to target Rac1 and inactivates it by increasing GTP hydrolysis [38]. Here we suggest a key role for miR-6715-3p in targeting Arha-gap12 and in this fashion, influencing Rac1 in the inner ear ( Figure 8).
Arhgap12, a GAP protein, is involved in controlling the time period in which Rac1 will remain active. In MDCK cells, Arhgap12 expression overlaps with Rac1 expression [38,39], probably due to the regulation of Arhgap12 on Rac1 at this time point. In the inner ear, besides the actin-rich areas, Rac1 is localized to hair cells and in particular, the kinocilium [37]. In accordance with our data, Arhgap12, at P0, with expression in the hair cells and kinocilium, might be regulating by Rac1 at this time, when the hair bundle is still developing and needs to be under close regulation. The kinocilium is necessary for proper development [40] and may have a structural role in guiding the architecture of the hair bundle and/or be involved in G-protein signaling in the kinocilium. If the latter is true, then Arhgap12 may regulate its GTPase activity. In our model of the miR-6715-3p circuit in the inner ear, we propose that the miRNA indirectly enhances the action of Rac1 and therefore may promote hair bundle development, actin reorganization, cell adhesion and inner ear morphogenesis.

Conclusions
The identification of inner ear-related miRNAs by RNAseq analysis demonstrates that the dataset is reliable not only for characterizing expression profiles of known miRNAs, but also for discovery of novel miRNAs in the inner ear. Further investigation of these miRNAs may shed light on their regulatory roles in various molecular pathways underlying the development of the embryonic inner ear.

RNA isolation
All procedures involving animals were approved and met the guidelines described in the National Institutes of Health Guide for the Use of Laboratory Animals and   approved by the Animal Care and Use Committee of Tel Aviv University (M-10-087). Cochlear and vestibular sensory epithelia were dissected from 12 inner ears of 6 P0 C57Bl/6 J mice. The cochlear sensory epithelia included the organ of Corti with some attached membranes, while the vestibular sensory epithelia included the saccule, utricle, lateral crista and anterior crista. The dissected sensory epithelia were placed immediately in Qiazol lysis buffer (Qiagen) on ice and stored at −80°C for several hours. The sensory epithelia was then thawed and total RNA was isolated, using the miRNeasy Mini Kit (Qiagen) without DNAse, which enables purification of total RNA of >18 nt. The isolated RNAs from the samples were precipitated in sodium acetate and glycogen in 80% ethanol (Ambion protocol), in order to further clean and concentrate the RNA samples. The precipitated RNA samples were dissolved in nuclease-free water and their quality was assessed by 1.2% agarose gel (90 V).

RNA-sequencing
Twelve cochlear epithelia were pooled and twelve vestibular epithelia were pooled, and each pool was used to create a library for MPS and sequenced. No biological repeats were performed. miRNA libraries for MPS were created using the TruSeq SmallRNA SamplePrep Kit (Illumina), with one modification. Instead of running the samples on acrylamide gel in order to isolate the miRNA band, the samples were loaded on a 4% agarose E-Gel and the purification was performed with the Qiagen Kit. The libraries were multiplexed and 2.5 pM of each library was sequenced at SR × 36 bp using the Illumina Genome Analyzer IIx in the same lane at the Functional Genomics Laboratory at Tel Aviv University.

miRNA sequence analysis
A total of 8,763,589 and 10,048,818 raw sequences of small RNAs were found in the cochlear and vestibular samples, respectively. Small RNA adapter sequences were clipped using fastq-mcf (ea-utils: http://code.google.com/ p/ea-utils/), discarding sequences shorter than 16 nucleotides and adapter dimers. Following adapter clipping, 8,340,741 and 9,441,679 reads remained. Differential expression of the resulting sequences was performed using miRNAkey. Specifically, sequence alignment was performed using Burrows-Wheeler Aligner (BWA) [41] against mouse mature miRNA sequences (downloaded from miRBase (Release 20: June 2013; http://www.mirbase.org/), while allowing one mismatch between read and reference. The read count for each read was normalized using the RPM method. Differential expression was tested using chi-squared proportion testing and corrected using the Bonferroni correction for multiple comparisons. Novel miRNA prediction was performed using miRDeep2 [42].

qRT-PCR
Cochlear sensory epithelia of E16, P0 and P8 C57BL/6 mice was dissected. Small RNAs were extracted using the miRNeasy Mini Kit (QIAGEN). Custom-made probes, designed by Applied Biosystems, were used to detect miR-NAs. miRNAs and U6B RNA (endogenous control) were reverse transcribed using the High Capacity cDNA Reverse Transcription Kit (Applied Biosystems). The qRT-PCR reaction was conducted using the FastStart Universal Probe Master (Roche) in the StepOne Plus qRT-PCR machine (Applied Biosystems). All miRNAs expression was normalized to the expression of U6B. At each age, sensory epithelia from 4-5 mice were pooled and five RT-PCR experiments were performed.

In situ hybridization
Inner ears of P0 C57BL/6 mice were dissected and fixed with 4% paraformaldehyde. Whole mount in situ hybridization analysis was performed using the Exiqon protocol. miRCURY LNA™ microRNA detection custommade probes, labeled with digoxygenin (DIG) (Exiqon), were used for the detection of novel miRNAs. Probes were hybridized with the tissue, 20-22°C below Tm of the probe. Probes were detected with the anti-DIG-AP (alkaline phosphatase conjugated) antibody (Roche), and the color reaction was developed using the NTB/BCIP (Sigma). The ears were frozen and cryosectioned to 10-18 μm sections and mounted. Images were taken using the Zeiss Aviovert200 M microscope. Three experiments were performed, and 3-5 ears were included in each experiment.

Luciferase assay
3' UTRs of chosen potential targets were cloned into the pGL3 luciferase reporter vector, downstream to the luciferase reporter gene, creating the 3'UTR-wt-pGL3 vector. The first three nucleotides of the miRNA binding site in the 3'UTR were mutated by site directed mutagenesis, creating the 3'UTR-mut-pGL3 vector. The miR-6715 miRNA expression vector was created by cloning the genomic polycistronic miR-6715-3p and miR-6715-5p into the miRvec expression vector (obtained as a gift from Reuven Agami) [43]. HEK293T cells were transiently transfected with either 3'UTR-wt-pGL3 or with 3'UTRmut-pGL3, for each of the tested targets, together with the corresponding miRNA expression vector and Renilla expressing vector. The luciferase reporter assay was performed 48 hours following transfection using the Dual-Luciferase® Reporter assay system (Promega). Three experiments were conducted, with duplicates.

Statistics
For qRT-PCR and luciferase assays, the Student's two-tailed t test P values of less than 0.05 were considered to be statistically significant and those of less than 0.005 were considered to be highly statistically significant. The data in the figures are presented by mean + SEM. Differential expression of miRNAs between paired samples was measured using a chi-squared statistic. P-values are calculated for the null hypothesis of no differential expression between the two samples. Final P-values were corrected using the Bonferroni correction for multiple hypotheses testing [18].