Developmental profiling of microRNAs in the human embryonic inner ear

Due to the extreme inaccessibility of fetal human inner ear tissue, defining of the microRNAs (miRNAs) that regulate development of the inner ear has relied on animal tissue. In the present study, we performed the first miRNA sequencing of otic precursors in human specimens. Using HTG miRNA Whole Transcriptome assays, we examined miRNA expression in the cochleovestibular ganglion (CVG), neural crest (NC), and otic vesicle (OV) from paraffin embedded (FFPE) human specimens in the Carnegie developmental stages 13–15. We found that in human embryonic tissues, there are different patterns of miRNA expression in the CVG, NC and OV. In particular, members of the miR-183 family (miR-96, miR-182, and miR-183) are differentially expressed in the CVG compared to NC and OV at Carnegie developmental stage 13. We further identified transcription factors that are differentially targeted in the CVG compared to the other tissues from stages 13–15, and we performed gene set enrichment analyses to determine differentially regulated pathways that are relevant to CVG development in humans. These findings not only provide insight into the mechanisms governing the development of the human inner ear, but also identify potential signaling pathways for promoting regeneration of the spiral ganglion and other components of the inner ear.


Introduction
MicroRNAs (miRNA) are a class of endogenously expressed small non-coding RNAs that function in RNA silencing and post-transcriptional regulation of gene expression.The human genome encodes more than 1000 miRNAs that may target as many as 60% of human protein-salaries for authors MS, BL, and BJK.MS had a role in Methodology and Investigation (performed HTG Edgeseq miRNA whole transcriptome assay), BL had a role in Software (wrote initial R code for our differential expression analysis), and BJK had a role in Project Administration.The specific roles of these authors are also articulated in the 'author contributions' section.

Competing interests:
We have read the journal's policy and the authors of the manuscript have the following competing interests: MS, BL, and BJK are employees of HTG Molecular Diagnostics, Inc., which has developed HTG Edgeseq miRNA whole transcriptome assay used in this study.HTG Molecular Diagnostics, Inc. does not alter our adherence to PLOS ONE policies on sharing data and materials.miRBase v20 database without the need for extracting RNA [19].Lysis of collected cells in each tube was performed by applying 55μL of HTG Lysis Buffer (HTG Molecular, Tucson, AZ) to the cap and incubating upside down for 30 minutes at room temperature.After 30 minutes the sample was spun down and 2.8μL of Proteinase K (HTG Molecular, Tucson, AZ, USA) was added, and the sample was then incubated for 180 minutes at 50˚C.From each prepared sample, 25μL were added per well to a 96-well sample plate.Human fetal Brain RNA was added to one well at 25 ng/well to serve as a process control.Samples were run on an HTG EdgeSeq Processor using the HTG EdgeSeq miRNA WTA (HTG Molecular, Tucson, AZ, USA).This assay revolves around nuclease protection, where a pre-selected miRNA Coronal section of a fetus with highlighted regions of interest (dashed black circles, CVG, NC, and OV).Further images are provided showing tissue before (C, E) and after (D, F) laser-captured microdissection of CVG (dashed white lines in D) and NC tissues (dashed white line in F).Other developmental regions also are highlighted.CVG: cochlear-vestibular ganglions, OV: otic vesicle, NC: neural crest; GG: geniculate ganglion, G-CVG: geniculatecochleovestibular ganglions, E: epithelium; and NT: neural tube.Scale bar = 100 μm (B) and 150 μm (C-F).
https://doi.org/10.1371/journal.pone.0191452.g001population is protected with proprietary protection probes, followed by degradation of all nonhybridized probes and non-targeted RNA by S1 nuclease.This step results in a 1:1 stoichiometric ratio of probes to targeted RNA.Following the processor step, samples were individually barcoded (using a 16-cycle PCR reaction to add adapters and molecular barcodes).Barcoded samples were individually purified using AMPure XP beads (Beckman Coulter, Brea, CA, USA) and quantitated using a KAPA Library Quantification kit (KAPA Biosystem, Wilmington, MA, USA).The library was sequenced on an Illumina MiSeq (Illumina, Inc., San Diego, CA) using a V3 150-cycle kit with two index reads.PhiX was spiked into the library at 5%; this spike-in control is standard for Illumina sequencing libraries.Data were returned from the sequencer in the form of demultiplexed FASTQ files, with one file per original well of the assay.The HTG EdgeSeq Parser (HTG Molecular, Tucson, AZ, USA) was used to align the FASTQ files to the probe list to collate the data.Data were provided as data tables of raw, quality control (QC) raw, counts per million (CPM), and median normalized counts.

Quality control
Baseline performance characteristics were established using Human Universal Reference RNA (uRNA) (Agilent genomics, Santa Clara, CA, USA) across all 96 wells on three sequencing runs of 96-well, with each plate processed on a different HTG EdgeSeq processor.Log 2 counts per million (log 2 (CPM)) standardization was used to transformed counts and adjusted for total reads within a sample [20].The RNA-sequence data consist of a matrix of read count r gi , for RNA samples i = 1 to n, and genes g = 1 to G. Ri, the total number of mapped reads for sample I, is defined as: The log 2 (CPM) values for each probe within a sample was defined as: Where y gi is normalized count for sample i and gene g.The counts are offset by 0.5 to avoid taking the log of zero, and total counts are offset by 1 to ensure the ratio is strictly less than 1 [20].
Baseline performance was established on the 96 technical replicates after normalization.Statistical process control methods were used to establish this expected performance where averaged ANT values (negative control) were calculated for each sample; a grand mean was calculated by taking the average of these averaged ANT values; the difference between each averaged ANT value and the grand mean, Δmean (averaged sample mean-grand mean) and standard deviation (SD) of the Δmean was calculated was for each sample.Acceptable Δmean values are those within ± 2SD average ANT.The graphical representation of the statistical process control, the quality control chart, for ANT using 96 miRNA technical replicates is shown in S1 Fig.

miRNA differential expression analysis and principal components analysis
A median normalization was used to transform and standardize data before data analyses [21].The primary statistical analysis was performed within the three time points and examined the differences in the pattern of expression among CVG, NC, and OV cell types.The R Bioconductor package (version 3.5) Linear Models for Microarray and RNA-Seq Data (limma) was used for the statistical analyses [22].Secondary statistical analysis examined differences in the pattern of expression over time within the CVG, NC, and OV cell types.The R Bioconductor package (version 3.5) time course was used for this analysis [23].The top 100 differentially expressed miRNA molecules across time points, ranked by significance (Hotelling's t-squared statistic), were determined for each of the three tissues.The maximal test statistic for the comparisons across three time points (stage 13 vs. 14, 14 vs. 15, and 13 vs. 15) was tabulated; differentially expressed miRNAs, indicating potential developmental relevance, were then identified by comparing these lists.Differential expression analyses included a correction for false discovery rate (FDR) was performed using the Benjamini-Hochberg Method [24].The R Bioconductor (version 3.5) was used for the analysis [25].To understand the primary determinant of miRNA expression and to examine sample clustering, Principal Component Analysis (PCA) was performed to identify the drivers of differences.Normalized miRNA expression values were assembled into a matrix with rows of different sample types and columns of miRNAs.

Over-connected gene analysis
Hypergeometric enrichment was used to identify genes significantly over-connected to miR-NAs detected in our differential expression datasets.Gene targets of miRNA were predicted by combining the TargetScan (version 7.1) nonconserved and conserved family databases, using all predicted targets for the human species [26].This provided targets for 1408 of the miRNA families that were used for over-connected gene analysis.A list was compiled of all miRNAgene links, allowing for multiple targets per miRNA.The hypergeometric distribution was used to compute p-values: where p is the probability that a gene would be targeted at least k times from a finite population under sampling without replacement if each gene had an equal probability of being targeted, N is the total number of miRNA-gene objects for all miRNA tested in the differential expression experiment, K is the number of times a given gene is represented in N, n is the number of miRNA-genes with parent miRNA that are considered differentially expressed (p < 0.05), and k is the number of miRNA-gene objects for a given gene represented by differentially expressed miRNA.This analysis was repeated for each miRNA differential expression dataset.The data were filtered for FDR corrected p-value < 0.10.Ensemble Gene Ontology (GO) annotations were obtained using the R Bioconductor (version 3.5) BiomaRt package [27].The complete list of over-connected genes was filtered to include only GO terms containing "transcription factor" in their annotation.To visualize these results, over-connected gene network graphics were created for several differential expression datasets using the R igraph package (version 1.0.1)[28].Significantly enriched miRNA-gene objects are represented with miRNA and genes as nodes.Node size reflects the number of connections made to each gene.Network graphics were created for all tissue comparisons with few enough over-connected transcription factors for the network plot to remain sparse enough to be useful.

Gene set enrichment analysis
The R Bioconductor mdgsa package (version 1.8.0) was used for pathway analysis of miRNA differential expression [29].R (version 3.3.3(Another Canoe), released on 03/06/2017) was additionally used for all subsequent statistical and computational analyses.The same TargetScan databases used for pathway analysis were used to predict miRNA-gene associations.Pathway-gene associations were taken from the Consensus Pathway Database (Consen-susPathDB, release 31), which combines 32 publically available interaction resources to provide information on protein interaction, signaling reactions, metabolic reactions, gene regulation, genetic interaction, drug-target interaction, and biochemical pathways [30].See further details in S1 File.
Integrated gene set analysis was performed following the method of Garcia-Garcia et al. [31], incorporating the paradigm shift proposed by Godard and Eyall to prevent the influence of knowledge bias [32].See details in S1 File.Following this approach, the pathway database (containing entries linking pathways to genes) was converted to an annotation database containing pathways linked to the miRNA targeting their genes using TargetScan predictions.Gene set analysis was then performed using mdgsa package, which fits a logistic regression model relating the probability of a gene belonging to the gene set with the value of the r statistic.The package outputs a log odds ratio for each interrogated gene set, along with raw and false discovery adjusted p-values [24].A positive log odds ratio indicates more targeting of a pathway by miRNA in case 1 compared to case 2 for each case 1 vs. case 2 comparison.
Alternatively, we identified significantly altered pathways using another type of gene set enrichment analysis.This analysis was conducted on the transformed gene expression from the miRNA analysis, which generated r statistics for each gene using the Garcia-Garcia method [31].The subsequent enrichment score is generated using R the Gene Set Variation Analysis (GSVA) package (version 1.24.1) on each gene set using ranking statistics similar to the Kolmogorov-Smirnov test described in Hanzelmann et al [33].Gene Ontology Pathway Database (Broad Institute, downloaded on 08/16/2017) was the gene set database, assayed in each sample for which it included a total of 5917 gene sets that are known to be involved in biological processes, cellular components and molecular functions [34,35].The package outputs an GSVA enrichment score for each interrogated gene set.The statistical significance is determined by p-values that are calculated using Baysian goodness of fit model using the R limma package (version 3.32.5).
To reduce functional redundancy within the Gene Ontology terms identified with gene set enrichment analysis with GSVA enrichment score, we next used REVIGO (http://revigo.irb.hr/), a web based program where sim Lin is Lin's sematic similarity measure [36], S (c 1 , c 2 ) is the set of common ancestors of terms c 1 , c 2, and p(c) is the probability of a term [37,38].The score is termed uniqueness in REVIGO.The previously calculated enrichment p-value for each term was not factored into the query since the aim was to reduce functionally overlapping or semantically redundant pathways; however, only significantly enriched terms were included in the query.Given that only one percent of randomly paired terms will be allowed at a uniqueness score of C = 0.53, a less stringent C = 0.7 was used to allow for a slightly semantically diverse list without significantly compromising the goal of reducing redundancy.

Computational prediction of the miR-183 family targeted genes
To investigate how the miR-183 family targets gene expression in human inner ear development, we performed computational target prediction for the miR-183 family (miRNA-96-5p, miR-182-5p, and miR-183-5p).As miRNA can have numerous targets, we first computationally predict the targets by using a miRNA prediction algorithm, TargetScanHuman (release 7.1: June 2016).TargetScan predicts biological targets of miRNAs by searching for the presence of conserved 8mer, 7mer, and 6mer sites that match the seed region of each miRNA [39,40].We then compared these results to the Deafness Variation Database (version 8.0), as this database provides a most comprehensive list of genetic variation in genes known to be associated with human deafness (http://deafnessvariationdatabase.org/download)[41,42].Relevant targets were thereby identified as the intersection of TargetScan predictions and genes associated with deafness as listed in Deafness Variation Database.

Differential miRNA expression among tissues for three developmental stages
Boxplots are shown in Fig 2 for expression of differentially expressed miRNA where CVG has the highest expression level at three different time points across the three tissue types, CVG, NC, and OV.These results show differentially expressed miRNA with p < 0.05 and false discovery rate (FDR) adjusted p < 0.25, 0.50, or 0.10 for tissue from Carnegie stages 13, 14, and 15, respectively.The smallest p-value out of the three possible tissue comparisons for each miRNA was considered.FDRs were controlled separately for each time point to minimize type II error while still maintaining a reasonable number of results, keeping in mind that a FDR controlled at 0.50 implies that 50% of differential expression results are expected to be false positives.Triplicate measurements were taken by collecting tissue from three different sections of one biological specimen.All three tissue samples within a time point are similarly from a single biological specimen.
Note that at the Carnegie developmental stage 13, members of the miR-183 family (miR-96, miR-182, and miR-183) were highly expressed in CVG and OV as compared with NC (for CVG vs NC comparison, p = 0.0009, 0.003, and 0.0045; with FDR correction, p = 0.27, 0.40, and 0.41, respectively).Also, note that expression of the entire miR-183 family was the highest in CVG tissue at stage 13 (see red circles on

Identification of over-connected transcription factors
To identify differential targeting of transcription factors, genes over-connected to miRNA differential expression datasets were identified by hypergeometric enrichment.The top 10 moretargeted and top 10 less-targeted transcription factors, ranked by significance and filtered for FDR adjusted at p < 0.10, are given for CVG vs. NC and CVG vs. OV comparisons for each

Gene set enrichment analysis
To understand tissue specific gene expression profiles that are important to human inner ear development, differentially expressed pathways were determined for each tissue comparison and time point based on miRNA differential expression datasets.The top 16 up-and downregulated pathways, ranked by significance, are included in Table 3 (CVG vs. NC) and 4 (CVG vs. OV), respectively.Note that a positive log odds ratio implies more miRNA targeting of a pathway, and therefore presumably inhibition of the pathway itself.In the CVG vs. NC comparison, the Dicer pathway, Insulin growth factor-1 (IGF-1) pathway, activated point mutants of FGFR2, and regulation of commissural axon pathfinding by SLIT and ROBO are likely relevant for CVG development (Table 3), because they have negative log odds ratios, and are thus less likely to be inhibited by miRNAs and more likely to be enriched in the ganglion as compared to the NC.In the CVG vs. OV comparison, regulation of commissural axon pathfinding by SLIT and ROBO is once again likely relevant for CVG development (Table 4).Complete differential pathway analyses data are included in the S5 Dataset.
Since each method for performing pathway analysis can have method dependent bias, we additionally performed a non-parametric analysis using the concept of gene set enrichment analysis [33] on the gene expression profile predicted from the miRNA expression pattern in CVG, OV and NC.From this analysis, we again identified functional gene sets defined in Gene Ontology terms that are positively or negatively enriched in each tissue type at different time points (Fig 7).analysis described in Supek et al., removing pathways that are functionally redundant after analyzing Gene Ontology classification [46] (Fig 8).Specific themes and processes were identified comparing Gene Ontology terms identified in CVG when comparing to NC at stage 13 (Fig 8 ) and at stage 14 (S6 Fig) .For example, six overarching biological processes were identified by examining enriched pathways at developmental stage 13 when comparing CVG to NC: regulation of lipoprotein metabolism, oligosaccharide biosynthesis, response to inactivity, hemidesome assembly, sarcoplasmic reticulum calcium ion transport, and replicative senescence.Note that replicative senescence includes pathways that are likely to relevant to CVG development such as auditory cell development and dorsal spinal cord development (red circled in Fig 8).A complete list of all enriched pathways is available in the S6 Dataset.

Computational prediction of the miR-183 family targeted genes
To further investigate how the miR-183 family targeted genes are expressed in the inner ear (CVG) development in humans, we computationally predicted miR-183 family targeted genes using TargetScanHuman (release 7.1).More than 1,000 genes were identified (complete datasets are included in S7 Dataset).The Deafness Variation Database (version 8.0) was then used to identify the targets that are associated with human hereditary hearing loss and deafness.Table 5 shows the list of predicted miR-183 family target genes that are relevant to human hereditary hearing loss and deafness.

Discussion/conclusions
Due to inaccessibility of human inner ear tissue, prior miRNA expression studies have been limited to animals, and in particular, mouse models [10,[47][48][49][50]. Further, the standard use in the United States and elsewhere of suction curettage (51), which generally renders embryonic structures unidentifiable, now makes it difficult even to obtain human tissue for study.However, the development of techniques for analyzing miRNA expression, along with the FFPE specimens provided by the Bruska group, enabled us to perform these novel analyses of human embryonic tissues.Note that due to the small sample size, our results should be cautiously assessed; bearing in mind that differential expression across time points can be confounded by biological noise, while differences among tissues could be specific to the sampled individual.
The differential expression of miR-183 family members (miR-96-5p, miR-182-5p, and miR-183-5p) in human CVG cell types at stage 13 correlates well with prior studies in mice.Weston et al. found that expression of the triad of miR-96, miR182, and miR-183 during development is relatively restricted to mouse inner ear compared to brain, heart, and whole embryo [3].Sacheli et al. noted the earliest expression of miR-183 and miR-182 in the mouse otic vesicle at embryonic day 9. (E9.5), expression of all three miRNAs in OV, CVG, and neural tube at E11.5, and limited expression was observed at E14.5 [51].By P0, this triad was strongly expressed in hair cells of the cochlea and vestibular system, as well as in the SGNs.We found that in human CVG, the miRNA-183 family was differentially expressed at stage 13 (E11), but only miR-182-5p was expressed at stage 14 (E11.5),and the miRNA family was not differentially expressed at stage 15 (E12).This may reflect differences between mouse and human development, but the small sample size in this study limits firm conclusions.Regardless, our findings suggest that the miR-183 family is important in human as well as murine inner ear development.
Twelve miRNAs out of the 100 most differentially expressed miRNAs across three time points were commonly expressed in all three tissues as presented in our Venn diagram (Fig 3).The majority of the common miRNAs (let-7 family, miRNA 548, 1255, 1273, 1285, 1301, and 1306) are involved in developmental timing, cell proliferation, cell-cycle regulation, and apoptosis [52,53].Also, miR-1254 is known to regulate the epithelial to mesenchymal transition (EMT) [54].EMT plays an important role in the delamination and escape through into the mesenchyme of both neural crest cells and placodal sensory neurons such as CVG [55].The commonality of expression of these twelve miRNAs among the CVG, OV, and NC cell types suggests that these miRNAs are involved in both inner ear and neural crest development.
The over-connected transcription factor analysis revealed a number of transcription factors that are known to be involved in CVG development in mice.For example, SKI was found to be upregulated in the CVG. vs. NC comparison, indicating more inhibition of SKI by miRNA in the CVG.SKI family transcription factors have been known to downregulate TGF-β signaling in mice thereby promoting neuronal induction including in the CVG [56].CARF (BDNF transcription, stage 13), EPHA5 (axonal genesis and synapse guidance, stage 13), LRP6 (Wnt/betacatenin signaling, stage 13), and BBS7 (Sonic hedgehog (SHH) signaling, stage 14) have been previously shown to be relevant in mouse CVG development [57].Other transcription factors identified by our analysis could represent important developmental effectors that have not been identified previously; our data thus has the potential to generate a series of follow up experiments to investigate development in the CVG, NC, and OV.
Our initial gene set enrichment analysis revealed that the Dicer, IGF-1, and FGFR2 pathways (stage 13) and regulation of commissural axon pathfinding by SLIT and ROBO (stage 14 and 15) were relevant ontogenetic pathways for CVG development.The Dicer pathway was the most statistically significant ontogenetic pathway at stage 13 in humans.A previous study   disruption of morphogenesis [4].More recently, a Dicer1 conditional KO line (Atoh1-cre; Dicer1 flox/flox and Foxg1-cre; Dicer1 flox/flox ) showed defects in proliferation in the prosensory domain of the cochlea including CVGs [7,8,58].Our observations coupled with the findings in mice suggest that Dicer pathway interactions with the miR-183 family play an essential role in the regulation of early mice inner ear development.Similarly, previous studies in mouse demonstrated that FGF2 signaling through FGFR2 and endogenous IGF-1 signaling are necessary for development of SGNs [59][60][61], consistent with our findings in the embryonic human inner ear.
In addition to identifying novel regulatory pathways derived from differential miRNA expression in specific embryonic inner ear tissues, we found that many of the enriched pathways occur at different stages of inner ear development.While comparing CVG to NC and OV, we observed that many more pathways were significantly enriched at stage 15, but that regulatory pathways of inner ear development, neurogenesis or axis patterning seemed to occur earlier at stages 13 and 14.Specifically, pathways involved in auditory receptor cell morphogenesis as well as many cell fate determination and pattering pathways were enriched in CVG at stages 13 and 14.By contrast, many of the pathways enriched in stage 15 are related to biosynthesis, metabolism and proliferation, suggesting that this is likely an embryonic stage allowing for growth after cell fate determination has already occurred.
Identification of miRNA-targeted genes is one of the most significant elements for understanding the function of a miRNA in the human inner ear development [47].Among our predicted miR-183 family targeted genes that are related to human inner ear (CVG) development (Table 5), it is worthwhile to mention TBX1 (T-box domain 1), as a previous study has confirmed tbx-1 as one of target genes for the conserved miR-183 family in a mouse model [62].Tbx1 was required for morphogenesis and growth of the murine otocyst [63].In humans, TBX1 is a critical gene in DiGeorge syndrome, with patients presenting phenotypes including inner ear abnormalities, sensorineural hearing loss [64][65][66], and vestibular loss [67,68].Therefore, we speculate that the transcriptional regulation of TBX1 in human may be influenced by the miR-183 family, providing a crucial function in developmental signaling pathways in the human inner ear.
It should be noted that this study is limited by small sample size; for each time point, only one biological specimen was obtained.It is, therefore, plausible that expression differences between time points reflected biological variations (i.e., noise) between individuals rather than tissue-specific variation over time.Differences between tissue types for a given time point could similarly be influenced by individual-specific variation.The small sample size also limits the statistical power to detect differential miRNA expression.While our findings are suitable for hypothesis generation, they should be critically considered on a case-by-case basis, bearing these limitations in mind.Additionally, though we would have liked to validate our findings using secondary methods such as qPCR, these experiments could not be performed due to sample constraints; all available tissue was used to obtain miRNA expression data with adequate quality control from HTG-seq.Many samples yielded only 500-1000 cells upon laser dissection (S1 Table ).Despite these limitations, many of our findings correlate well with prior findings in the mouse which suggests that they are truly informative about inner ear development in humans.In addition to providing insight into the mechanisms governing the development of the human inner ear, our findings also identify potential signaling pathways for promoting regeneration of human spiral ganglions and other components of the inner ear using pluripotent stem cells in the future.
Fig 2).The trend was not seen at Stage 14, where the only member of the miR-183 family that was upregulated was miR-182 (p = 0.004; with FDR correction p = 0.27) (indicated with a green circle in Fig 2).At stage 15, no members of the miR-183 family were observed.Boxplots for expression of the top-20 differentially expressed over the three time points for CVG, NC, and OV are shown in S2 Fig.The top-100 differentially expressed miRNAs in CVG, OV, and NC, and the overlap across tissues, are summarized in Fig 3A as a Venn diagram.A total of twelve miRNAs were found to be commonly regulated over time in these three tissues; let-7d, let-7e, miR-1249, miR-1254, miR-1255, miR-1273, miR-1285, miR-1301, miR-1306, miR-548, and miR-8078.Differential expression analyses at three developmental time points on these twelves miRNAs are shown in Fig 3B.The complete set of raw data obtained from HTG Edgeseq miRNA WTA, Quality control files, Brain RNA correction, and normalized data are shown in the S1 Dataset.Complete spreadsheets of differential miRNA expression for each tissue comparison and time point are available in the S2 and S3 Datasets, respectively.Principal component analysis (PCA) was performed to visualize potential clustering among samples by either tissue or time point.Principal components 1 and 2 are plotted against one another in Fig 4.Although there was no tight clustering among samples by tissue or time point, PC2 is moderately associated with tissue type (Fig 4A, dashed circles).Components are alternatively visualized as the means among technical replicates for a given tissue type and time point with x and y standard error bars (Fig 4B).

Fig 2 .
Fig 2. Differential expression of miRNA across tissues.Differentially expressed miRNA between three tissue types (CVG, NC, and OV) with p < 0.05 and false discovery adjusted p < 0.25, 0.50, or 0.10 (based on the smallest p-value out of the three possible tissue comparisons) are plotted for tissue from Carnegie stages 13, 14, and 15, respectively.MicroRNA from the miR-183 family is circled in red (stage 13) or green (stage 14).https://doi.org/10.1371/journal.pone.0191452.g002

Fig 3 .
Fig 3. Differentially expressed miRNA across time points.(A) A Venn diagram is presented showing the top 100 miRNAs differentially expressed among the three time points (stages 13, 14, and 15) for each tissue type (CVG, NC, and OV).(B) Normalized counts are plotted for 12 commonly expressed miRNAs that are differentially expressed across time for all three tissues.https://doi.org/10.1371/journal.pone.0191452.g003 Fig 7 shows two signatures of the differentially activated Gene Ontology pathways by GSVA analysis in CVG vs. NC comparison (Fig 7A) and CVG vs. OV (Fig 7B) in heatmap format.Using this method, we found that different numbers of pathways were identified as significantly enriched when comparing different tissue types at different stages of development (S5 Fig), with stage 15 having highest number of enriched pathways (209 pathways) when comparing CVG to NC or OV.To generate a signature, any pathway that was significantly enriched either positively or negatively at p-<0.05 in at least one time point was included in the heatmap.Note that there are many GO term pathways directly relevant to CVG development (Fig 7C and 7D).These include Wnt activated receptor activity, ganglion development, axon extension, positive regulation of neuron migration, and integrin mediated cell adhesion (higher GSVA score in Stage 15) at different time points (Fig 7D) while most common category of Gene Ontology terms for all comparisons are pathways involved in general biological processes (shown in blue bars in S5 Fig).To better focus on the biological processes that are enriched in different tissue types, we additionally performed the REVIGO

Fig 4 .
Fig 4. Principal component analysis.(A) Principle components 1 and 2 are plotted for the 27 collected datasets.Loose clustering by tissue type is indicated with dashed colored circles.(B) As an alternative visualization, triplicate measures for each individual tissue/time point set are averaged and plotted with x-and y-standard error bars.https://doi.org/10.1371/journal.pone.0191452.g004

Fig 5 .
Fig 5. Transcription factor targeting network plots, CVG vs NC comparison.Network plots showing differentially targeted transcription factors between CVG and NC tissues.Green nodes represent differentially expressed miRNAs, while orange and purple nodes represent significantly (false discovery adjusted p < 0.10) targeted transcription factors more targeted in either CVG (purple) or NC (orange) tissue.Transcription factor node size is proportional to the number of differentially expressed miRNAs targeting the gene.https://doi.org/10.1371/journal.pone.0191452.g005

Fig 6 .
Fig 6.Transcription factor targeting network plots, CVG vs OV comparison.Network plots showing differentially targeted transcription factors between CVG and OV tissues.Green nodes represent differentially expressed miRNAs, while orange and purple nodes represent significantly (false discovery adjusted p < 0.10) targeted transcription factors more targeted in either CVG (purple) or OV (orange) tissue.Transcription factor node size is proportional to the number of differentially expressed miRNAs targeting the gene.Only results for stages 13 and 15 are shown, as there are too many significant factors to effectively visualize at stage 14. https://doi.org/10.1371/journal.pone.0191452.g006

Fig 7 .
Fig 7. Differentially activated gene ontological pathways by GSVA analysis.Heatmap comparison of CVG with NC (A), and of CVG with OV (B).Each column represents the sample indicated at the top (three time points), each row represents an identified GO Term pathway (see the complete set in the S6 Dataset).The expression levels of GSVA scores are depicted according to the color scale (middle right).Red or green indicate expression levels above or below the median, respectively.The magnitude of deviation from the median is represented by color saturation.Heatmaps representing the GSVA score of GO term pathways relevant to CVG development that are significantly modulated in CVG. vs. NC comparison (C) and CVG vs. OV (D).https://doi.org/10.1371/journal.pone.0191452.g007

Fig 8 .
Fig 8. REVIGO treemap of relevance similarity analysis on enriched pathways comparing CVG vs. NC (stage 13).REVIGO treemap summarizing Gene Ontology biological process categories over-represented in CVG cells compared to NC cells at stage 13.All terms are included with a FDR adjusted p-value cutoff at 0.05 from the enrichment analysis.The relevance similarity C-score (uniqueness) cut-off is chosen at 0.7.The size of each rectangle is proportional to the uniqueness for that category.Red circles indicate relevant CVG-development pathways.https://doi.org/10.1371/journal.pone.0191452.g008 Table.Detailed description of the three FFPE samples.C-Stage: the Carnegie developmental stage; CR length: Crown-rump length; NC: neural crest; CVG: cochlear-vestibular ganglions; OV: otic vesicle.(DOCX) S1 Fig.The quality control chart.SD: standard deviation.The difference between averaged ANT (control) median normalized counts and the grand mean of the averaged ANT values (Δmean) was computed for each sample and plotted to demonstrate statistical process control for HTG measurement techniques.Samples were deemed acceptable if they fell within ±2SD.(PDF) S2 Fig. Boxplots for expression of the top 20 differentially expressed miRNAs over the three time points for CVG, NC, and OV.Plots of normalized counts versus time for the top 20 differentially expressed miRNAs for each tissue ranked by significance.Error bar = ±2SD, S#1: sample 1; S#2: sample 2; S#3: sample 3. (PDF) Table 5. Predicted target genes for the MiR-183 family using TagetScanHuman and a Deafness Variation Database (version 8.0).miRNA Predicted target genes miR-183 family MITF, TECTB, EDNRB, CD164, OSBPL2, MET, CLIC5, GRHL2, COL4A4, TBX1, RDX, S1PR2, CCDC50, PNPT1, PEX6, SLC26A4, OTOG, COL9A1, COL4A6, PTPRQ https://doi.org/10.1371/journal.pone.0191452.t005S3 Fig. Heatmap of differentially activated gene ontological pathways by GSVA analysis (CVG vs. NC).(TIFF) S4 Fig. Heatmap of differentially activated gene ontological pathways by GSVA analysis (CVG vs. OV).(TIFF) S5 Fig. Enriched pathways by GSVA analysis.Comparisons among three tissue types at three time points.Blue bar: biological function; orange bar: molecular function; gray bar: cellular component; and yellow bar: others.The number of pathway is shown inside or next to each bar.S: Stage.(PDF) S6 Fig. REVIGO treemap of relevance similarity analysis on enriched pathways comparing CVG vs. NC (stage 14).REVIGO treemap summarizing Gene Ontology biological process categories over-represented in CVG cells compared to NC cells at stage 14.All terms are included with a FDR adjusted p-value cutoff at 0.05 from the enrichment analysis.The relevance similarity C-score (uniqueness) cut-off is chosen at 0.7.The size of each rectangle is proportional to the uniqueness for that category.Red circles indicate relevant neuronal development pathways.(PDF) S1 Dataset.Raw expression count data.Raw expression data for each sample as provided to us by HTG-seq.Raw counts for each sample are converted to counts per million (CPM), then median normalized in the "Normalized" tab.Spiked in brain RNA is used as a control, expecting high correlation between groups.Samples are matched to their normalized data in the "Groupings" tab, then averaged prior to computing fold changes between tissues in the "Fold Change" tab.(XLSX) S2 Dataset.Differential expression data across three tissue types.Differential expression data based analysis with the R Bioconductor limma package.The spreadsheet includes averaged normalized miRNA expression values for each tissue, fold changes for all tissue comparisons, p-values, and false discovery adjusted p-values from limma.Each time point is included as a separate tab.(XLSX) S3 Dataset.Differential expression data across three time points.Differential expression data analyzed using the R Bioconductor timecourse package.The spreadsheet includes averaged normalized miRNA expression values for each time point and the maximum Hotelling T values from the three possible time point comparisons.Results are ordered by rank.Each time point is included as a separate tab.(XLSX) S4 Dataset.Differentially targeted genes.Complete differential gene targeting analysis for each tissue comparison and time point (separate tabs), as analyzed by hypergeometric enrichment.Complete results are included, as well as filtered results for each comparison at the false discovery adjusted p < 0.10 level.Genes that contain "transcription factor" in their Gene Ontology are highlighted in light green.(XLSX) S5 Dataset.Differentially regulated pathways.Complete gene set analysis results for each tissue comparison and time point (separate tabs), analyzed using the R Bioconductor mdgsa