Identification of microRNAs associated with human fragile X syndrome using next-generation sequencing

Fragile X syndrome (FXS) is caused by a mutation in the FMR1 gene which can lead to a loss or shortage of the FMR1 protein. This protein interacts with specific miRNAs and can cause a range of neurological disorders. Therefore, miRNAs could act as a novel class of biomarkers for common CNS diseases. This study aimed to test this theory by exploring the expression profiles of various miRNAs in Iranian using deep sequencing-based technologies and validating the miRNAs affecting the expression of the FMR1 gene. Blood samples were taken from 15 patients with FXS (9 males, 6 females) and 12 controls. 25 miRNAs were differentially expressed in individuals with FXS compared to controls. Levels of 9 miRNAs were found to be significantly changed (3 upregulated and 6 downregulated). In Patients, the levels of hsa-miR-532-5p, hsa-miR-652-3p and hsa-miR-4797-3p were significantly upregulated while levels of hsa-miR-191-5p, hsa-miR-181-5p, hsa-miR-26a-5p, hsa-miR-30e-5p, hsa-miR-186-5p, and hsa-miR-4797-5p exhibited significant downregulation; and these dysregulations were confirmed by RT‐qPCR. This study presents among the first evidence of altered miRNA expression in blood samples from patients with FXS, which could be used for diagnostic, prognostic, and treatment purposes. Larger studies are required to confirm these preliminary results.

www.nature.com/scientificreports/ (Eppendorf, Hamburg, Germany) using the AmplideX PCR/CE FMR1 Kit (Asuragen, Austin, TX) and Ampli-deX PCR/CEFMR1reagents (cat. no. 49402), following the manufacturer's instructions. PCR products were then separated in a 3310XL capillary electrophoresis system [Applied Biosystems Genetic Analyzer, Foster City, CA, USA] based on conditions described in the kit manual. GeneMapper software (ID-X Software Version 1.0), Applied Biosystems, was used to analyze and convert the separated PCR products into CGG repeat length using.
In accordance with the current ACMG Guidelines, the CGG triplets that repeated 45-54 and 55-200 times were considered as intermediate and premutation respectively. Those with more than 200 repeats were defined as full mutation. Samples with both premutation and full mutation alleles were identified as full mutation mosaics.
Total RNA quantity and quality control. All isolated RNA samples were eluted in RNase-free water and RNA concentrations were determined with Nanodrop Spectrophotometer (Thermo Scientific, Wilmington, DE, USA). RNA quantity was evaluated by calculating absorbance at λ = 260 nm, and the quality was assessed by a ratio of λ = 260/280 nm being close to 2.0-2.3. The RNA concentration of each sample was more than 50 ng. The integrity of the RNA, as a key feature that affects the performance of sequencing and RT-qPCR, was assessed via two methods: First, by running extracted RNA through 1% agarose gel and then staining with ethidium bromide to observe the 28S ribosomal RNA band at 4.5 kb, and the 18S rRNA band at 1.9 kb. Second, by using the Agilent 4200 TapeStation System (Agilent Technologies, Santa Clara, CA, G2991A) to assess the electrophoretic profile of the 18S and 28S RNA and generating an RNA Integrity Number (RIN). All RNA samples revealed RIN values of greater than eight and the miRNA extractions were stored at − 80 °C until processing.
Library preparation for next-generation sequencing. The TruSeq Small RNA Library Preparation Kit (Illumina, San Diego, USA) was used for generating miRNA sequencing libraries directly from total RNA as per the manufacturer's protocol for this kit (TruSeq Small RNA Library Prep Reference Guide 15004197 v02). Briefly, after ligation of the 5′ and 3′ RNA adapters using T4 RNA ligase, reverse transcription was performed to generate cDNA; cDNA libraries were subsequently amplified by PCR. The products were then purified. After acrylamide gel purification, eight libraries were pooled in an equimolar amount to create one lane and validation was compatible with multiplexed sequencing. Finally, the libraries were checked and normalized according to protocol.
miRNA profiling through next-generation sequencing (miRNA-seq). The miRNA cDNA libraries were sequenced on an Illumina MiniSeq platform in the Pars gene company, Shiraz, Iran. With this platform, DNA fragments of the libraries go through clonal amplification by bridge PCR followed by sequencing using a reversible terminator. It consists of sequencing by synthesis (SBS) technology using only two-channel (i.e., red and green) which needs only two images for determination of all four base calls reducing the number of cycles, cost, and time required for data processing and yet delivering high accuracy and quality 23 . Within every cycle of sequencing, for each cluster, base calls were created by inbuilt real-time analysis software and raw data were stored in the format of individual base call files (*.bcl). The BCL files were converted to standard FASTQ file formats (a text-based sequencing data file format that stores both raw sequence data and quality scores) for downstream analysis. FASTQ data acquisition in available repositories, Sequence Read Archive (SRA: Accession number: PRJNA777620) in NCBI was obtained and is in processing to transfer to GEO.
In the next step, the output data were streamed into Illumina's BaseSpace Sequence Hub for cloud-based data management and analysis. FASTQ files were cleaned by adapter removal using CutAdapt 1.6. The Phred numerical quality scoring system was used as a base call quality filter 24,25 ; reads with Q scores of < 33 (checked both before and after adapter trimming) were removed. After removal, the adaptors and filtering out low-quality sequences, the processed reads were aligned against the miRBase database and human genome hg19 by using version 2.2.3 of Short Read Mapping Package (SHRiMP). The number of reads mapped to miRNAs in each sample provided as SHRiMP log table and trimming-mapping plot are illustrated in supplementary 1, 2 respectively. According to the SHRiMP log table, the range of Reads is 519045 to 2456946 and the percentage of mapped reads for all miRNAs is 2.79-83.67% as minimum and maximum. Mapped files were then sorted and indexed as binary format (BAM) files.

Relative quantification of miRNA by reverse-transcription PCR analyses. Quantitative reverse
transcription-polymerase chain reaction (qRT-PCR) is a well-established method for miRNA profiling with the highest sensitivity and accuracy and with the widest dynamic range 26 . We used stem-loop RT-PCR miRNA assay for quantification of miRNA expression levels 27 using a commercially available qRT-PCR miRNA detection kit and primer sets (Zist Fanavari Pishgam Company, Tehran, Iran). To normalize the expression data, a commonly-used endogenous reference gene, U6 small nuclear RNA (U6-snRNA), was adopted 28 . Real-time PCR data with the use of a LightCycler 96 instrument (Roche Applied Science, Indianapolis, IN, USA), were analyzed by 2 −ΔΔCT method using the following equation: ΔΔCT = ((Ct miRNA of concern − Ct U6-snRNA) in patients) − ((Ct miRNA of concern − Ct U6-snRNA) in controls) 29 . These relative gene expression data analyses were double-checked by using the relative expression software (REST-2009) tool 30 .
The MicroRNAs gene expression was compared between sexes and number of repeats (genotypes) and was analyzed using independent sample t and ANOVA tests (Table 5).
To evaluate the degree of similarity between RT-qPCR and miRNA-seq results, we compared miRNAs' expression data generated through both approaches. Downregulation was defined when 2 −ΔΔCT was less than zero and upregulation presented as more than zero results of 2 −ΔΔCT . www.nature.com/scientificreports/ Differentially expressed miRNA analysis of deep sequencing data. To find differentially expressed miRNAs between case and controls, we have to compare expression levels of nearly 1800 miRNAs that have read count numbers greater than zero for at least one of the two groups. This leads to 1800 statistical tests. When multiple hypotheses are tested simultaneously, an adjustment for the multiplicity of tests is often necessary to restrict the total number of false discoveries, and usage of such methods has become standard in genomics 31 . The most conservative method for adjusting p values is Bonferroni correction 31,32 . There exist other adjustment methods which increase the power of the statistical test in the cost of decreasing control over type I error.
Our study using Bonferroni correction leads to over a hundred candidate miRNAs but we just want to choose a handful number of them for further investigations, so statistical power is not an issue here.
A reasonable way of choosing miRNAs for further exploration would be choosing miRNAs with the lowest p values. However, because the lowest p values vary significantly by adding or removing some samples, they may be too sensitive to noise. Also, their differences are too small and their ordering may not be robust too. So we used a slightly different approach other than simply picking miRNAs with the lowest p values. In this study, to determine differentially expressed miRNAs, we conducted differential expression analysis using four different software tools, namely, DESeq2 version1.26.0 (https:// bioco nduct or. org/ packa ges/ relea se/ bioc/ html/ DESeq2. html), edgeR version 3.28.1 (https:// bioco nduct or. org/ packa ges/ relea se/ bioc/ html/ edgeR. html), limma-trend and limma-voom version 3.42.2 (https:// bioco nduct or. org/ packa ges/ relea se/ bioc/ html/ limma. html). We listed 30 miRNAs with the lowest p values in each method and computed intersections of these 30 candidates. Interestingly, the agreement among the methods in calling differentially expressed miRNAs was high, despite their differences in assumptions and algorithms and also in the order and magnitude of their resulting p values. The overall intersection consists of 25 miRNAs which showed noteworthy agreement between different methods (Table2).

Result CGG allele sizing and clinical characteristics of studied subjects. The clinical and demographic
findings of the studied patients are listed in Table 1. The specific genotype category results are characterized as full mutation (n = 10), premutation (n = 3), and full mutation mosaic (n = 2). In this study, all (100%) available full mutation and full mutation mosaic patients with FXS presented with ASD, and one-third (33%) had ADHD.
Differentially expressed deep sequencing of miRNA. As explained in "Differentially expressed miRNA analysis of deep sequencing data" section, we first picked 25 miRNA in (Table2), utilizing four differential expression analysis tools DESeq2, edgeR, limma-trend, and limma-voom. To apply further investigations on the miRNAs we needed to filter this list to obtain only 9 miRNAs. The number of chosen miRNAs depends on what type of further exploration we are going to apply and its cost.
These 25 miRNAs cannot be sorted based on p values because each of the DE analysis tools used, imposes a different p value ordering on this list. From the perspective of fold changes, the absolute value of fold changes for all the resulted miRNAs was over 2 with a majority of them over 3, which is not surprising, because fold change is already somehow taken into account when computing p values. Hence we did not use p value or fold change as the parameter for this round of filtering.
We sorted the list of these 25 miRNAs based on the overall mean expression level for all cases and controls and chose the first nine highly expressed miRNAs. The intuition behind this selection was that we expected www.nature.com/scientificreports/ these miRNAs to present more powerful signals and lower noises as opposed to miRNAs with a lower level of expressions and we expect their p values to be more robust. miR-125, though was not in 30 highly dysregulated candidate miRNAs in this study, was also analyzed due to its importance in urine in the Putkonen et al. 19 study and the existence of a small number of NGS studies about FXS. As seen in Table 3, all members of the miR-125 family were found to be downregulated in patients with FXS.
Proposed miR reconfirmation by RT-PCR. The schematic heat map shows log expression-related changes in miRNA transcriptome based on three groups of full mutations, permutation, full mutation mosaic (Fig. 1).
Regarding the difference between the expression of miRs in the three groups of premutation, full mutation, and mosaic, only miR-532-5p has a significantly different distribution in the three groups. Using ANOVA, the most overexpression was found in full-mutation and then premutation and mosaic respectively (p value = 0.0423, Df = 2, F-value = 4.165).  www.nature.com/scientificreports/ There was no difference between men and women in genes expression of miRs (p vlaue > 0.05). Some miRNAs including miR-191-5p, miR-26a-5p, and miR-181-5p were downregulated, whereas miR-652-3p was upregulated in all six females. These miRNAs seemed to be gender-dependent. Similarly, age dependency was seen in those aged over 55 years with miR-191-5p, miR-26a-5p, and miR-181-5p. Among These miRs, expression of miR-532-5p gene was significantly more in patient under 55 years than over 55 years (p value = 0.03416, df = 4.5305, t = 2.9979).
In addition, no significant correlation has existed between any miRNAs and ADHD (p value > 0.05). The presence of seizure disorder was documented in three patients that all showed downregulation of miR-26a-5p and miR-186-5p. Autism spectrum disorder was significantly present in full mutation patients with downregulated miR-181-5p (p value < 0.05). In search in targetscan7.7 database, (http:// www. targe tscan. org/ vert_ 72/) for prediction of miRNA target, no target was found for miR-181-5p.
Finally, the end result of relative expression with 2 −ΔΔCT calculation are summarized in Table 5 and

Discussion
The current study identified FXS-specific changes in miRNAs among Iranian blood samples as preliminary evidence. We identified twenty-five differentially expressed miRNAs sequenced in the blood of individuals with FXS compared to the controls, and we found minor downregulation of miR-125a-5p.
MicroRNAs regulate mRNAs at the post-transcriptional level and therefore affect protein translation 12 . Changed miRNA expression patterns epigenetically affect almost every aspect of CNS function (i.e., in neurogenesis, synaptogenesis, and neuronal migration) and its development [33][34][35] . For instance, miR-532 is reliably expressed in the human brain, localized as distinct granules in distal axons and growth cones, and proposed to play a role in axon growth and guidance 36 . ZFHX3 gene, among 5 target genes of the hsa-miR-532-5p from the MiRTarBase microRNA Targets dataset, encodes a transcription factor that regulates neuronal differentiation 37 .
FXS is the first neurodevelopmental disease found to be linked to the dysfunction of the miRNA pathway 17 . The in vivo evidence of miRNA involvement in FXS pathogenesis was first provided in a study of the zebrafish model by identifying and isolating numerous miRNAs, including miR-FMR1-27 and miR-FMR1-42 in this model 38 .
Subsequent studies in the FMR1 KO mouse models found that disruption of the regulating of miR-125a, miR-125b, and miR-132 causes early neural development and synaptic physiology 18 and that there is an interaction between miR-34b, miR-340, and miR-148a with the Met 3′ UTR of the FMR1 gene 39 . Moreover, by isolating mesenchymal stem cells from peripheral blood and differentiating these cells into neuronal cells, Fazeli et al. 40 recently analyzed the expression of miR-510 by the qPCR method. The authors reported an enhanced expression of miR-510, located on chromosome X in the 27.3Xq region, flanking to a fragile X site, in the female carriers of FMR1 full mutation. www.nature.com/scientificreports/ In a most recent study, Putkonen et al. 19 showed upregulation of miR-125a in urine from children with FXS. The investigators did not examine differential miRNA expression changes in FXS blood samples or the correlations of miR-125a levels in urine with those of in the cell-free circulation (i.e., in serum and plasma) and other body fluids 19 . In our study, we found a minor downregulation of miR-125a-5p in the blood of individuals with FXS. One preliminary hypothesis for this finding is that due to urinary secretion of miR-125a-5p, its blood level expression decreased similar to what happens for the blood-urine balance of electrolytes. In line with our results in the Alvarez-Mora MI study on FXTAS patients using deep sequencing, the authors observed a slight but not significant reduction of miR-125a-5p in the blood of FXTAS patients 41 . Recently, Frye et al. published an article on 2021 that found hsa-miR-125a-5p_R-1 was significantly down-regulated in the ASD group, which is in parallel to our current study 42 . They introduced the role of miR-181 in immunomodulation in ASD in lymphoblastoid cell lines derived from 10 individuals with ASD. They also found that out of 267 detected miRs in the ASD vs. Sibling group, 14 other miRNAs were dysregulated 42 . Furthermore, in a microarray analysis by Seno et al. the up-regulation of miR-125b as well as miR-196a, miR-650, miR-338-3p, in 20 cases of ASD were reported 43 .
Moreover, in concordance with our results, the association between autistic traits and X-linked SNPs in the gene family linked with FXS, is likely to be owing to a disruption in the recognition between has-miR-181 and the corresponding seed match sequences in these genes 45 ; miR-181d and FMRP cooperatively regulate the axon elongation process 22,46 . An altered expression pattern for miR-181 and miR-191 in hippocampal neuron development has been shown to occur 47 .
Even though larger studies are needed to confirm our results and investigate the effect of other miRNAs, the changes in miRNAs seen among our patients provide evidence that these miRNAs could have roles in developmental processes, nervous system homeostasis and the function of nerve cells in those with FXS.
In one recent study, 13 miRNAs were differentially expressed in maternal plasma samples from pregnant women with fetal Down syndrome versus healthy control subjects; among the others, hsa-miR-191 was upregulated and hsa-miR30e downregulated 48 . In another study, miR-26b-5p, miR-185-5p, and miR-191-5p were identified as potential biomarkers for ADHD in peripheral blood mononuclear cells 49 . Altered expression of miR-26a and miR-26b have been shown in peripheral blood of major depressive patients during antidepressant therapy, in Alzheimer's disease, and Parkinson's disease 50 . Finally, hsa-miR-532-5p and hsa-miR-652-3p are upregulated in schizophrenia 51,52 . www.nature.com/scientificreports/ Our finding demonstrates significant involvement of hsa-miR-30e-5p in FXS, which was found to be the most significantly upregulated miRNA in patients with FXS compared with controls. miR-30 family plays a major regulating role in the tissue and organ development and the pathogenesis of various clinical diseases 53 . Several studies have shown that hsa-miR-30e-5p among other miRNAs might be associated with the onset and progression of Parkinson's disease and schizophrenia [54][55][56]  Our results also showed deregulated hsa-miR-191-5p. Although, as far as we are aware, no evidence for hsa-miR-191-5p contribution to FXS has been reported so far, alterations in the expression level of hsa-miR-191-5p have earlier been found in patients with neuropsychiatric disorders sharing genetic overlap with FXS, including ASD, ADHD, schizophrenia, bipolar disorder, and major depressive disorder 49,50,61,62 .
It is noteworthy that despite a large number of miRNAs associated with FXS and its related disorders that have been identified in multiple expression studies, only a few miRNAs are common between various studies. This discrepancy can be explained in part by the polygenic and complex nature of neuropsychiatric disorders 63 .
The expression profiles of the miRNAs in our study confirm some existing findings but conflict with others. For example, our result regarding the expression level of has-miR-30e in patients with FXS is consistent with Sun et al. 's 60 report that miR-30e was upregulated in PBMCs from patients with schizophrenia. In contrast, Perkins et al. 55 have found that miR-30e is downregulated in the prefrontal cortex of subjects with schizophrenia compared with healthy subjects. The exact reason for these conflicting results remains to be determined but it may be because of differences in screening standards (i.e., patients' ethnicity, geographical region, and screening criteria), techniques used for miRNA detection and profiling, and experimental design. Furthermore, as previously mentioned by Alvarez-Mora et al. 41 , it is documented that the expressions of miRNAs are tissue-specific and/or temporally regulated which may partially explain the differences seen between the findings of different studies.
This study is based on a small sample size and is an initial discovery on the path of diagnostic, prognostic, and treatment purposes in FXS. Several types of research to identify biomarkers of relevance to clinical trials of targeted therapeutics in FXS recommended by experts in the field 64 . Some promising efforts also help address an ongoing issue of placebo using parent-based outcome measures in the clinical trials of targeted therapeutics in FXS [65][66][67] .
Several limitations to our study should be addressed. First, the sample size is relatively small. Despite this, it is among the first pieces of evidence of altered miRNA expressions in blood samples from patients with FXS. Second, although patients consisted of individuals from all across the country, they were only Iranian in origin. Third, Validation was performed on the same dataset, due to the small number of patients who were willing to take blood. Therefore, it does not provide the exact effectiveness of the miRNAs as biomarkers. It should be validated on an independent set of patients and donor specimens. Forth we examined miRNA expression changes in non-neuronal cells since neuronal tissue is not easily accessible. However, it has been shown that miRNA expression changes in the peripheral circulation are highly correlated with those of neuronal tissue from patients with various neuropsychiatric disorders 33 .

Conclusions
Our study is among the first to present the characterization of the miRNA's expression profiles in blood samples of patients with FXS using deep sequencing-based technologies. Out of 25 picked miRNAs, levels of nine miRNAs were found to be changed (i.e., 3 upregulated and 6 downregulated) in FXS blood. Altered peripheral miRNA levels have been identified in numerous neuropsychiatric disorders, including FXTAS, ASD, ADHD, Down syndrome, depression, and schizophrenia. Our results provide a new perspective for the role of miRNA profiling in the pathophysiology of FXS, but larger studies are required to confirm these preliminary results and explore the influence of the other dysregulated miRNAs. If confirmed, it could open the possibility of using miRNAs as novel non-invasive FXS biomarkers or broad-spectrum therapeutic agents.