MicroRNA-Related Polymorphism and Their Association with Fibromyalgia

MicroRNAs are tissue-specific expressed short RNAs that serve post-transcriptional gene regulation. A specific microRNA can bind to mRNAs of different genes and thereby suppress their protein production. In the context of the complex phenotype of fibromyalgia, we used the Axiom miRNA Target Site Genotyping Array to search genome-wide for DNA variations in microRNA genes, their regulatory regions, and in the 3’UTR of protein-coding genes. To identify disease-relevant DNA polymorphisms, a cohort of 176 female fibromyalgia patients was studied in comparison to a cohort of 162 healthy women. The association between 48,329 markers and fibromyalgia was investigated using logistic regression adjusted for population stratification. Results show that 29 markers had p-values < 1 × 10−3, and the strongest association was observed for rs758459 (p-value of 0.0001), located in the Neurogenin 1 gene which is targeted by hsa-miR-130a-3p. Furthermore, variant rs2295963 is predicted to affect binding of hsa-miR-1-3p. Both microRNAs were previously reported to be differentially expressed in fibromyalgia patients. Despite its limited statistical power, this study reports two microRNA-related polymorphisms which may play a functional role in the pathogenesis of fibromyalgia. For a better understanding of the disease pattern, further functional analyses on the biological significance of microRNAs and microRNA-related polymorphisms are required.


Introduction
Fibromyalgia (FM) is a chronic pain disorder characterized by widespread musculoskeletal pain and tenderness at specific sites [1]. Patients report symptoms at the joint muscles and tendons on all four quadrants of the body, with increased pain under load. Additional symptoms include general weakness, chronic fatigue, sleep disturbances, cognitive dysfunction, impaired concentration, and reduced mental and physical capacity [2,3]. Patients with FM have reduced pain thresholds and impaired neuronal processing and modulation of pain [4][5][6][7]. The prevalence of fibromyalgia is estimated at around 2-4% in the general population (85-90% are female) and may be caused by a combination of genetic and environmental factors [8][9][10]. Certain environmental factors, such as physical or emotional trauma, may trigger the development of fibromyalgia in individuals who are genetically predisposed to the condition [11][12][13].
MicroRNAs (miRNAs; miR), first described in 1993 [14], are small non-coding single-stranded RNA molecules that play an essential role in gene regulation at the The data collection took place between May 2019 and July 2020. The study was approved by the ethics review board of the Medical Faculty, Ruhr University Bochum (6594-BR). All participants gave written informed consent to participate in the study.

DNA Extraction and Quality Controls, and Genotyping Array
DNA was extracted from mouthwash samples using the salting-out method as described by Miller et al. [38]. DNA concentration and purity was photometrically measured (Synergy 2, Biotech/Agilent; Santa Clara, CA, USA) and integrity estimated by agarose gel electrophoresis. A total of 376 samples (197 patients with FM, 179 HC) had the required yield, purity, and integrity, were diluted to a concentration of 15 ng/µL in a volume of 80 µL, and were sent to our partner laboratory Cologne Center for Genomics (CGC). DNA of all subjects was simultaneously analyzed for a total of 237,858 DNA variations distributed in the genome using Axiom miRNA Target Site Genotyping Arrays (Applied Biosystems/Thermo Fisher Scientific; Waltham, MA, USA). This genome-wide array focuses on genotyping SNPs and InDels located in miRNA promoters, miRNA seed sequences, in the 3'UTR of protein-coding genes and in the DNA of miRNA-processing proteins. In a final step, Axiom miRNA Target Site Genotyping Arrays were processed using the GeneTitan ® Multi-Channel Instrument, according to the manufacturer's protocol (Thermo Fischer Scientific, Waltham, MA, USA).

Biostatistical Analysis
The generated array data were first pre-processed using Thermo Fisher Axiom Analysis Suite software. Further quality control was performed using Plink [39]. Here, standard quality control steps were performed using the following exclusion criteria for individuals: >0.02 missingness, heterozygosity rate > |0.20| or sex-mismatch. SNPs were excluded if they deviated from Hardy-Weinberg Equilibrium (HWE) with a p-value of <10 −6 , had minor allele frequency <0.01 or missing data >0.02. Relatedness and population stratification were determined using a SNP set filtered for SNP quality (HWE p > 0.02, MAF > 0.20, missingness = 0) and LD adjustment (r 2 = 0.1). In case of cryptic relatedness (pi-hat > 0.20), one subject was randomly removed. To control for population structure, subjects who differed by more than 4.5 standard deviations in the first 20 genetic principal components (PCs) were excluded. As a result, data of 48,329 SNPs were available from 170 fibromyalgia patients and 162 healthy controls.
SNPs related to FM were identified by applying a genome-wide association study (GWAS) approach. In particular, SNPs were tested for association using the logistic regression model incorporated in Plink. Principal components 1-5 from quality control were thereby included as covariates. Bonferroni-corrected genome-wide significance was calculated as p < 1.03 × 10 −6 (0.05/48,329 SNPs). Though according to the strategy implemented by Wilkins and colleagues [40], SNPs with p-values < 1 × 10 −3 were also considered for analysis of miRNA-related functional consequences.
In silico analysis of predicted miRNA-binding alterations due to DNA-polymorphisms were performed by mining the miRNASNP-v3 database [41]. miRNASNP-v3 makes use of TargetScan v7.2 and miRmap v1.1 to determine miRNA target binding sites and to predict SNP-induced gain or loss of miRNA/target pairs [41][42][43]. Accordingly, we examined all of the p < 1 × 10 −3 SNPs for inducing gain or loss of miRNA/target pairs. Subsequently, we checked whether any of the SNP-induced gains or losses of miRNAs were previously related to FM. A total of 39 FM-related miRNAs were identified from manually scanning 26 publications reporting on miRNA expression in FM. The literature search was performed on PubMed based on the query "miRNA AND fibromyalgia" on the 20 February 2023. Furthermore, we made use of miRDB [44,45]; (target prediction score ≥ 95) to predict gene targets of 37 of the 39 FM-related miRNAs. Targets of hsa-miR-133a and hsa-miR-146a could not be predicted due to lack of specification whether the 3p or 5p miRNA was investigated in the original publication. Ultimately, we checked for overlaps between the predicted gene targets and genes harboring the p < 1 × 10 −3 SNPs. Lists of miRNAs identified from the literature (see Table S1), their predicted targets (see Table S3) as well predicted SNP-induced gain or loss of miRNA/target pairs (see Table S2) are available in the Supplementary Materials.

Genome-Wide DNA-Polymorphism Associations
The test for association between 48,329 miRNA-related SNPs and fibromyalgia revealed 29 SNPs with p-values below our significance threshold of p < 1 × 10 −3 , with rs758459 displaying the smallest p-value of 0.0001 (see Figure 1, Table 2). makes use of TargetScan v7.2 and miRmap v1.1 to determine miRNA target binding sites and to predict SNP-induced gain or loss of miRNA/target pairs [41][42][43]. Accordingly, we examined all of the p < 1 × 10 −3 SNPs for inducing gain or loss of miRNA/target pairs. Subsequently, we checked whether any of the SNP-induced gains or losses of miRNAs were previously related to FM. A total of 39 FM-related miRNAs were identified from manually scanning 26 publications reporting on miRNA expression in FM. The literature search was performed on PubMed based on the query "miRNA AND fibromyalgia" on the 20th of February 2023. Furthermore, we made use of miRDB [44,45]; (target prediction score ≥ 95) to predict gene targets of 37 of the 39 FM-related miRNAs. Targets of hsa-miR-133a and hsa-miR-146a could not be predicted due to lack of specification whether the 3p or 5p miRNA was investigated in the original publication. Ultimately, we checked for overlaps between the predicted gene targets and genes harboring the p < 1 × 10 −3 SNPs. Lists of miRNAs identified from the literature (see Table S1), their predicted targets (see Table S3) as well predicted SNP-induced gain or loss of miRNA/target pairs (see Table S2) are available in the Supplementary Materials.

Genome-Wide DNA-Polymorphism Associations
The test for association between 48,329 miRNA-related SNPs and fibromyalgia revealed 29 SNPs with p-values below our significance threshold of p < 1 × 10 −3 , with rs758459 displaying the smallest p-value of 0.0001 (see Figure 1, Table 2).    Note. Additional information is received from the Thermo Fisher Axiom miRNA Target Site Genotyping Array annotation or from dbSNP 156. SNPs are sorted by their p-value. * According to the Thermo Fisher annotation, rs2295963 resides in SPATA6L. However, dbSNP 156 additionally labels rs2295963 as a PLPP6 3-Prime UTR Variant.

Predicted miRNA/Target Pair Gain or Loss
Mining of the miRNASNP-v3 database revealed that 16 out of the 29 SNPs were predicted to alter binding of altogether 171 unique miRNAs (see Table S2). Comparison with miRNAs identified from the FM literature exhibited an overlap with hsa-miR-1-3p (MI0000437) whose binding was predicted to be affected by rs2295963 at the Phospholipid Phosphatase 6 gene (PLPP6). Precisely, presence of rs2295963 predicted loss of the hsa-miR-1-3p/PLPP6 pair. Additionally, rs2295963 was predicted to induce loss of 7 and gain of 4 additional miRNA/PLPP6 pairs (see Table 3, Figure 2).

Overlap of Genes and Predicted miRNA Targets
A number of 1245 unique genes were predicted as functional targets from 37 literatureextracted miRNAs by miRDB (see Table S3). Comparison with genes of p < 1 × 10 −3 SNPs revealed an overlap with the Neurogenin 1 gene (NEUROG1), which was predicted to be targeted by hsa-miR-130a-3p and contained the SNP with the lowest p-value (rs758459; see Table 2).

Overlap of Genes and Predicted miRNA Targets
A number of 1245 unique genes were predicted as functional targets from 37 literature-extracted miRNAs by miRDB (see Table S3). Comparison with genes of p < 1 × 10 −3 SNPs revealed an overlap with the Neurogenin 1 gene (NEUROG1), which was predicted to be targeted by hsa-miR-130a-3p and contained the SNP with the lowest p-value (rs758459; see Table 2).

Discussion
FM is a complex chronic pain syndrome characterized by individually varying clinical characteristics and a substantial genetic component [46]. Former GWAS of chronic pain and FM identified multiple SNPs residing in genes that affect central nervous system functioning [47,48] and calcium regulation [47][48][49]. Moreover, an increasing number of miR-NAs has been associated with FM pathogenesis (for an overview, see [50]). However, there is still lack of evidence regarding SNPs altering miRNA production and binding in FM. Although SNPs in pre-miRNA loci yield the potential to disrupt production of mature miRNAs, and SNPs in miRNA target genes may alter binding of miRNAs (reviewed by

Discussion
FM is a complex chronic pain syndrome characterized by individually varying clinical characteristics and a substantial genetic component [46]. Former GWAS of chronic pain and FM identified multiple SNPs residing in genes that affect central nervous system functioning [47,48] and calcium regulation [47][48][49]. Moreover, an increasing number of miRNAs has been associated with FM pathogenesis (for an overview, see [50]). However, there is still lack of evidence regarding SNPs altering miRNA production and binding in FM. Although SNPs in pre-miRNA loci yield the potential to disrupt production of mature miRNAs, and SNPs in miRNA target genes may alter binding of miRNAs (reviewed by [51]), none have been focused on in FM research.
We have therefore made use of the Axiom miRNA Target Site Genotyping Array to investigate miRNA-related SNPs on a genome-wide level in 170 FM patients and 162 HC. More than 80% of the SNPs that can be genotyped using this array are not found on other genotyping arrays. Accordingly, utilizing this miRNA SNP array was particularly interesting for our research question. To compensate the reduced statistical power that results from our small sample size, we conducted in silico analyses and evaluated our results based on the literature on miRNA-expression in FM. According to our hypothesis, SNPs with p-values < 1 × 10 −3 should possibly have pleiotropic effects (e.g., by affecting miRNA binding sites), influencing molecular signaling pathways and protein expression. In fact, one of the most promising SNPs (rs2295963) identified by our GWAS was predicted to alter binding of hsa-miR-1-3p, a miRNA known to be differentially expressed in FM patients [27]. Moreover, the in silico analysis proposed 11 more miRNAs whose binding may be altered due to rs2295963, suggesting presence of pleiotropic effects. Larger case control studies could replicate this association and give more precise estimates of the effect sizes. Since rs2295963 is located in the 3'UTR for phospholipid phosphatase 6 (PLPP6), a gene that has not yet been described in connection with FM, it would be the logical next step to investigate the effects of this SNP on mRNA stability and PLPP6 protein expression. Luciferase assays and addition of miRNA mimics and antagomirs to cell cultures could provide further insights into the possibilities of the regulatory influence of this SNP. In a further step, it would be interesting to explore which pathways show involvement of PLPP6 and to what extent and under which conditions these pathways contribute to the pathogenesis of FM. In addition, studies of free circulating miRNAs during different phases of pain intensity could provide new insights on how specific miRNAs behave in acute and chronic pain as it was already shown for cell-free DNA under stressful life conditions [52]. These studies could be complemented by drug administration to better understand pharmacological influences on the release of miRNAs in FM pathogenesis.
Further comparisons of predicted miRNA targets and genes of SNPs passing the p-value < 1 × 10 −3 threshold detected an overlap with NEUROG1 harboring rs758459. Even if NEUROG1 initially appears to be a promising candidate here, since it assumes important functions as a basic helix-loop-helix (bHLH) transcription factor in the context of neurogenesis, a miRNA-mediated function of rs758459 is probably ruled out here due to its localization in the promoter region of NEUROG1. However, a function in the context of NEUROG1 promoter activity cannot be ruled out and should therefore be examined more closely, e.g., by promoter assays as described elsewhere [53].
Nonetheless, our results are clearly limited by statistical power as well as the number of assessed SNPs. The total sample size of 332 subjects did not allow for discovery of genomewide significant polymorphisms. Moreover, odds ratios of the p < 1 × 10 −3 SNPs were high due to the small sample size. Additionally, the 48,329 SNPs investigated in this study represent only a fraction of all miRNA-related SNPs [41]. Future studies should therefore not only consider increasing the sample size but also the number of polymorphisms possibly altering miRNA binding and production, allowing for a comprehensive exploration of the miRNA-related genetic variation associated with FM. SNPs in coding regions of FM-related miRNAs are of special interest here, as they may provide an explanation for the already noted aberrant expression patterns of those miRNAs. Moreover, it is important to note that expression of miRNAs is dynamic and tissue-specific [54]. Genome-wide expression studies have shown that miRNAs are expressed differently in postmortem brain tissue depending on psychiatric phenotypes. In vitro studies confirmed that these miRNAs can modulate the protein expression of their specific target genes [55,56]. This is accomplished by epigenetic signatures and can be modulated by external factors [57,58]. Accordingly, results of miRNA-based studies could be further integrated into multi-omics studies to gain a better understanding of pathological altered gene regulation, starting from transcriptional regulation to protein expression [59][60][61][62].

Conclusions
To conclude, we conducted a GWAS of miRNA-related SNPs in a small FM cohort and identified two novel SNPs based on the respective miRNA expression literature. While rs2295963 was predicted to alter binding of hsa-miR-1-3p as well as binding of 11 additional miRNAs, rs758459 is located within a gene that is targeted by hsa-miR-130a-3p. Even though findings are preliminary, they reveal two promising polymorphisms for functional follow-up studies and a scalable, literature-driven approach to identify further miRNArelated SNPs that may contribute to the complex genetic architecture of FM.

Institutional Review Board Statement:
The study was conducted in accordance with the Declaration of Helsinki and was approved by the ethics review board of the Medical Faculty, Ruhr University Bochum (6594-BR).

Informed Consent Statement:
Written informed consent has been obtained from the patient(s) to publish this paper.

Data Availability Statement:
The data presented in this study are available on request. Individual genetic data cannot be made available due to privacy restrictions. The R scripts used are available on the publicly accessible OSF page of the corresponding author.