In silico identification and in vitro expression analysis of breast cancer-related m⁶A-SNPs

ABSTRACT

Research on m⁶A-associated SNPs (m⁶A-SNPs) has emerged recently due to their possible critical roles in many key biological processes. In this sense, several investigations have identified m⁶A-SNPs in different diseases. In order to gain a more complete understanding of the role that m⁶A-SNPs can play in breast cancer, we performed an in silico analysis to identify the m⁶A-SNPs associated with breast cancer and to evaluate their possible effects. For this purpose, we downloaded SNPs related to breast cancer and a list of m⁶A-SNPs from public databases in order to identify which ones appear in both. Subsequently, we assessed the identified m⁶A-SNPs in silico by expression quantitative trait loci (eQTL) analysis and differential gene expression analysis. We genotyped the m⁶A-SNPs found in the in silico analysis in 35 patients with breast cancer, and we carried out a gene expression analysis experimentally on those that showed differences. Our results identified 981 m⁶A-SNPs related to breast cancer. Four m⁶A-SNPs showed an eQTL effect and only three were in genes that presented an altered gene expression. When the three m⁶A-SNPs were evaluated in the tissue sample of our breast cancer patients, only the m⁶A-SNP rs76563149 located in ZNF354A gene presented differences in allele frequencies and a low gene expression in breast cancer tissues, especially in luminal B HER2+ subtype. Future investigations of these m⁶A-SNPs should expand the study in different ethnic groups and increase the sample sizes to test their association with breast cancer and elucidate their molecular function.

KEYWORDS:

• Breast cancer
• epigenetics
• m⁶A methylation
• m⁶A-SNP
• SNP

Introduction

Breast cancer (BC) is a heterogeneous disease, both biologically and clinically. Worldwide, BC is the most commonly diagnosed cancer among women [1]. Its incidence, mortality, and survival rates vary considerably between different geographical areas, which could be due to many factors such as population structure, lifestyle, genetic factors, and environment [2]. Genetic factors play an important role in the aetiology of BC. Investigations carried out by genome-wide association studies (GWAS) have successfully identified many genetic markers like Single Nucleotide Polymorphisms (SNPs) that are associated with BC risk [3]. Moreover, epigenetic changes have recently emerged as markers since they might be involved in cancer initiation and progression. Consequently, they might be useful for early detection, prognosis, and targeted therapy of cancer [4].

Epitranscriptomics is still an unexplored area that, like epigenomics, may be of great interest to respond to the relationship between environmental factors and cancer, since it involves modifications in the RNAs that modulate gene expression in response to external factors. In this sense, N⁶-methyladenosine (m⁶A) is the most common post-transcriptional modification in eukaryotic messenger RNA (mRNA) [5]. m⁶A is a reversible modification that occurs throughout the entire transcriptome into a consensus motif known as DRA*CH (D = G, A, or U; R = G, or A; A* = A methylatable, and H = A, C, or U) [6,7] mainly by the action of methyltransferases METTL3, METTL14, and WTAP (‘writers’), and can be reversed by the enzymes FTO and ALKBH5 (‘erasers’). In addition, multiple proteins, known as ‘readers,’ specifically bind to the methylated m⁶A position. The m⁶A plays a critical role in controlling mRNA structure, stability, splicing, export, localization, and translation [5, 8–10].

Recent works indicate that SNPs located within or near m⁶A motifs may affect m⁶A methylation and subsequent related biological processes by altering the RNA sequence of the target site or key flanking nucleotides [11,12]. These genetic variations are known as m⁶A-associated SNPs (m⁶A-SNPs). In this context, these m⁶A-SNPs could have an eQTL (expression Quantitative Trait Loci) effect. The cis-eQTL effect on the expression of local genes may result from multiple causes, such as altering the level of m⁶A methylation, the binding of transcription factors and proteins, the mRNA stability, nuclear export, transcription elongation, splicing, or modifying the secondary structure of mRNA [11, 13–15] and, consequently, influence the BC development and progression.

To date, several investigations have identified the role of m⁶A-SNPs in rheumatoid arthritis [16], osteoporosis [17,18], coronary artery disease [19], periodontitis [13], development of sepsis [14], adiposity [20], lipid levels (cholesterol and triglycerides) [21], Parkinson’s disease [22], inflammatory bowel disease [23], colorectal cancer [24], type 2 diabetes [25], and thyroid cancer [26]. Recently, Xuan et al. (2021) [27] have found eight m⁶A-SNPs that could have an effect on BC. In order to have a more complete understanding of the role those m⁶A-SNPs could play in BC, the present study aimed to carry out an independent parallel in silico analysis using publicly available data to identify new m⁶A-SNPs associated with this disease and evaluate their possible effects on gene expression.

Materials and methods

Determination of m⁶A-SNPs related to breast cancer

In order to determine the effect of the m⁶A-SNPs on BC, first, we had to identify the SNPs related to this disease. We searched for these SNPs through the NCI Genomic Data Commons (GDC, https://portal.gdc.cancer.gov/) [28], the GWAS Catalogue (https://www.ebi.ac.uk/gwas/) [29], and the SNPedia (https://www.snpedia.com/) [30]. Additionally, we also searched the SNPs present in the microRNAs (miRNAs) related to BC in the following databases: the miRdSNP (http://mirdsnp.ccr.buffalo.edu/index.php) [31], the miRNASNP (http://bioinfo.life.hust.edu.cn/miRNASNP/#!/) [32], the MiRNA SNP Disease Database (MSDD, http://bio-bigdata.hrbmu.edu.cn/msdd/home.jsp) [33], and the SomamiR (http://compbio.uthsc.edu/SomamiR/) [34]. Subsequently, a list of m⁶A-SNPs and the associated data were downloaded from the m6Avar database (http://m6avar.renlab.org/help.html) [12]. This database presents human m⁶A-SNPs with high, medium, and low confidence levels. Finally, we compared the list of the SNPs associated with BC and the m⁶A-SNPs with the purpose of identifying which ones appear in both.

Data on allele and genotype frequencies for the identified m⁶A-SNPs were obtained from Ensembl Genome Browser (https://www.ensembl.org/index.html) [35].

Expression quantitative trait loci analysis for the m⁶A-SNPs related to breast cancer

To find out if the identified m⁶A-SNPs change the expression of local gene, we carried out a search of cis-eQTL analyses using public data from the GTEx Portal (https://gtexportal.org/home/) and the HaploReg browser (http://archive.broadinstitute.org/mammals/haploreg/haploreg.php) [36]. For that, we introduced the m⁶A-SNPs identified in the databases and observed if they presented an cis-eQTL effect in the breast tissue when we used the data from GTEx Portal and in different human tissues and cells when we used HaploReg browser. Moreover, we also searched for potential functions in gene expression regulation, such as altering protein binding or the regulation motif in HaploReg browser.

Differential expression analysis

For the m⁶A-SNPs that showed cis-eQTL effect, we investigated if their genes presented a differential expression in BC using the expression profile data available in GEO database (http://www.ncbi.nlm.nih.gov/geo) [37]. For this purpose, we downloaded two datasets that contained the normalized gene expression signals for breast tissue (normal and tumour). The dataset GSE5764 contained data on gene expression levels of 10 individuals with BC (five with invasive ductal and five with invasive lobular breast carcinoma) and 20 corresponding to normal tissue (10 to normal ductal and 10 to normal lobular tissue) [38]. Dataset GSE21422 contained data from five healthy tissue samples, nine mammary ductal carcinomas in situ, and five invasive ductal carcinomas [39]. We analysed the differences in gene expression levels using the Mann–Whitney U-test through IBM SPSS Statistics 25.0 (SPSS Inc., IBM Corporation, New York, USA). We considered differential expression genes if they had p-values <0.05 and if they had been detected in at least one of the two studies. In the case of the m⁶A-SNPs in the miRNA that presented cis-eQTL effect, data about the number of samples, gene expression mean, and p-values were obtained from the ENCORI for RNA Interactomes database (http://starbase.sysu.edu.cn/panMirDiffExp.php) [40].

Validation of the m⁶A-SNPs in patients with breast cancer

Patients

In order to validate the m⁶A-SNPs that showed a cis-eQTL effect and differences in gene expression after in silico analysis, we analysed breast tissues from BC patients. Samples and associated clinicopathologic data were provided by the Basque Biobank after approval by the Basque Clinical Research Ethics Committee (PI2018044). Frozen breast tumour and corresponding adjacent normal-appearing tissues (located at least 3 cm away from the site at which the tumour was sampled) were collected from 35 patients. All patients presented invasive ductal BC, a mean age of 59.77 ± 8.63, and 68.57% of them were menopausal. The tumour samples were classified into five molecular subtypes: seven samples belonged to luminal A, seven to luminal B HER2-, seven to triple-negative breast cancer (TNBC), eight samples to luminal B HER2+, and six to HER2+ subtype. Tissue histology, as well as molecular subtype classification, was confirmed by two experienced pathologists.

Total RNA was extracted from the tissue samples using TRIzol® Reagent (Invitrogen, Carlsbad, USA) and PureLink® RNA Mini Kit (Invitrogen) following the manufacturer’s instructions. The cDNA synthesis was performed through SuperScript® VILO™ Synthesis kit (Invitrogen) following the manufacturer’s instructions.

As a control group of gene frequencies, data of Iberian populations in Spain (IBS) for the m⁶A-SNPs were downloaded from the 1000 Genomes database (https://www.internationalgenome.org/home) [41]. Only data from females were selected for this analysis.

Genotype determination

The genotypes for each m⁶A-SNP identified in the in silico analysis were determined for our patient samples by sequencing. The primers for each m⁶A-SNPs were designed using the PerlPrimer v1.1.21 program [42]. All primer pairs were verified using NCBI’s Primer-BLAST software (https://www.ncbi.nlm.nih.gov/tools/primer-blast/) and synthesized by Metabion (Planegg, Germany). The sequences of the primers, as well as the size of the amplified, and the annealing temperature can be seen in Table S1. We performed Real-time PCR (qRT-PCR) amplification on a C1000 Thermal Cycler coupled to a CFX96 Real-Time PCR Detection System (Bio-Rad, California, USA). Reactions were carried out in a final volume of 5 μL and contained: 0.50 μL of each forward and reverse primers (1 μM), 2.50 μL of SsoFast EvaGreen Supermix (Bio-Rad), 1.00 μL of cDNA template at 2.50 ng/μL, and 0.50 μL of milliQ water (Merck Millipore, Germany). The qRT-PCR conditions were as follows: an initial denaturalization at 98°C for 2 min, followed by 35 cycles at 98°C for 10 s and the annealing temperature for each set of primers for 30 s. Subsequently, the amplicons were purified with ExoSAP-IT^TM Express (Applied Biosystems, Waltham, USA). We performed sequencing reactions with BigDye® Terminator^TM v1.1 Cycle Sequencing Kit (Applied Biosystems) and we purified sequencing products with BigDye® XTerminator^TM Purification Kit (Applied Biosystems), following the manufacturer’s instructions. We performed capillary electrophoresis on an ABI PRISM 3130 Genetic Analyzer (Applied Biosystems). We analysed sequencing reaction products with the Sequencing Analysis software v.5.2 (Applied Biosystems) and we determined the genotype of the patients for each m⁶A-SNPs by ChromasPro software v.2.6.5 (Technelysium Pty Ltd, Brisbane, Australia).

Gene expression analysis

We carried out a gene expression analysis for the m⁶A-SNPs that presented significant differences in terms of genotype and/or allele frequencies. We conducted qRT-PCR of gene expression described in the previous section. All qRT-PCR analyses were performed in triplicate for each sample. The relative gene expression levels were estimated through the 2^−ΔΔCt method [43] using as reference genes GAPDH, β-ACTIN, and GUSB. We choose these reference genes since they are the most commonly used in this type of studies. However, recently, some authors reported variability between different tissues and conditions for these genes [44], which would make them not optimal as reference genes. Therefore, as other authors did to get over this drawback, we used an average of their Ct values [45,46]. Table S1 shows information about the sequence, annealing temperature, amplicon size, and efficiency of the reference genes primers.

Statistical analysis

Allele and genotype frequencies were determined for each m⁶A-SNP. Differences between the group of patients with BC and the control group from the IBS population were analysed using Pearson’s chi-square (χ²) test. On the other hand, the differences in gene expression between the tumour breast tissue of the patients and their corresponding adjacent normal breast tissues were analysed using the Mann–Whitney U-test. The variation between different tumour subtypes was determined using the Kruskal–Wallis H test. We used IBM SPSS Statistics 25.0 software (SPSS Inc., IBM Corporation, NY, USA) for all tests. We took values of p < 0.05 as statistically significant except in the Kruskal–Wallis H test where a Bonferroni correction was applied (with a statistically significant cut-off value of p < 0.005).

Results

m⁶A-SNPs related to breast cancer

In total, we identified 981 BC-related m⁶A-SNPs after comparing 69,713 SNPs related to BC (65,331 from GDC, 1,775 from GWAS Catalogue, 71 from SNPedia, 1,864 from miRdSNP, 42 from miRNASNP, 43 from MSDD, and 587 from SomamiR) with 301,875 m⁶A-SNPs in RNA (301,529 m⁶A-SNPs present in mRNA and 346 in miRNA) (Figure 1). From these 981 m⁶A-SNPs, there were 26 that belonged to the high confidence category, 147 belonged to the medium confidence category, and 808 to the low confidence category. These different confidence categories relate to the identification method for the m⁶A-SNPs [12]. PA-m⁶A-seq (photo-crosslinking-assisted m⁶A-seq) and miCLIP (m⁶A individual-nucleotide-resolution crosslinking and immunoprecipitation) methods provide high confidence data. On the other hand, MeRIP-seq (m⁶A-specific methylated RNA immunoprecipitation with next-generation sequencing) provides medium-confidence. Finally, a genome-wide prediction based on the Random Forest algorithm provides low confidence.

Figure 1. Flow chart of study design and the main results. IBS: Iberian Population in Spain from 1000 Genomes Project.

Of the 981 m⁶A-SNPs, all but one were in mRNA. The other one was in miRNA. It is well known that most m⁶A methylated sites are in the 5’ untranslated region (UTR), in long exons, near the stop codon, in the 3’-UTR, and non-coding RNAs (ncRNA) [8, 47–49]. In our study, we found that most of the m⁶A-SNPs, i.e., SNPs located within or near m⁶A motifs, are distributed in the coding sequence (90.62%), 6.63% are in the 3’-UTR, 2.65% in the 5’-UTR, and only 0.10% in the miRNA.

Cis-eQTL analysis

To further explore the potential effect of the 981 m⁶A-SNPs in BC, we investigated whether they present cis-eQTL effect. Considering the data from the GTEx Portal for BC samples, none of these m⁶A-SNPs presented a cis-eQTL effect. However, when we searched in the HaploReg database, four m⁶A-SNPs (rs1801270, rs9379084, rs76563149, and rs11614913) showed a cis-eQTL effect in different cells or tissues (Figure 1).

The m⁶A-SNP rs1801270 (C/A) presents cis-eQTL effect in fibroblasts [50]. This m⁶A-SNP is a missense variant localized in the Cyclin-Dependent Kinase Inhibitor 1A gene (CDKN1A, gene ID: 1026) on chromosome 6 (Table 1). Variant A of this m⁶A-SNP is found in an incomplete DRACH motif (AGA*CG, where the A* is the methylated adenosine). Moreover, this m⁶A-SNP is near one complete DRACH motif and one incomplete one (at 16 nucleotides upstream and at five nucleotides downstream in the RNA, respectively). The frequency of the A allele varies depending on the population, between 0.067 and 0.488 for European (EUR) and East Asian (EAS) populations, respectively (Table 2). On the other hand, according to the ENCODE transcription factor ChIP-seq datasets [51] obtained in HaploReg database, this m⁶A-SNP could also be involved in gene expression regulation by altering the regulatory motif, i.e., a binding site of transcription factors, and changing the binding site of proteins POL2 and ZBTB33.

Table 1. List of identified m⁶A-SNPs associated with BC that presented cis-eQTL effect.

m^6A-SNP	CHR^a	Position^b	Gene	Mutation type	Strand	Proteins bound	Motifs changed	m^6A ID
rs1801270	6	36,684,194	CDKN1A	Missense	+	POL2 and ZBTB33	Ascl2, ELF1 HEN1, and Sin3AK-20	m6A_ID_71975
rs9379084	6	7,231,610	RREB1	Missense	+	No	VDR	m6A_ID_72448
rs76563149	5	178,730,572	ZNF354A	5’-UTR^c	-	POL2, SIN3AK20, POL24H8, TAF1, E2F1, HAE2F1, TR4, CCNT2, HEY1, HMGN3, and ZNF263	EBF and Ets	m6A_ID_187781
rs11614913	12	53,991,815	MIR196A2	miRNA^d	+	No	HMG-IY	m6A_ID_358107

^aCHR: chromosome; ^b Assembly: GRCh38.p13; ^c 5’-UTR: 5’-untranslated region; ^d miRNA: microRNA

Table 2. Alleles and genotypes frequencies for the four m⁶A-SNPs in different populations. AFR: African; AMR: American; EAS: East Asian; EUR: European; and SAS: South Asian.

		Frequency of alleles and genotypes
m^6A-SNP	Alleles/Genotypes	AFR	AMR	EAS	EUR	SAS
rs1801270	AA + AC	0.526	0.421	0.716	0.127	0.262
	CC	0.474	0.579	0.284	0.873	0.738
	A	0.309	0.249	0.488	0.067	0.145
	C	0.691	0.751	0.512	0.933	0.885
rs9379084	AA + AG	0.012	0.216	0.284	0.247	0.202
	GG	0.988	0.784	0.716	0.753	0.798
	A	0.006	0.111	0.159	0.132	0.106
	G	0.994	0.889	0.841	0.868	0.894
rs76563149^a	TT + GT	0.156	0.412	0.357	0.508	0.593
	GG	0.844	0.588	0.643	0.492	0.407
	T	0.081	0.235	0.201	0.298	0.368
	G	0.919	0.765	0.799	0.702	0.632
rs11614913	CC + TC	0.977	0.838	0.714	0.817	0.930
	TT	0.023	0.162	0.286	0.183	0.070
	C	0.860	0.611	0.458	0.589	0.742
	T	0.140	0.389	0.542	0.411	0.258

^aThe variants for this m⁶A-SNP in the mRNA are A or C.

The other m⁶A-SNP is rs9379084 (G/A) that presented a cis-eQTL effect in whole blood. This m⁶A-SNP is a missense variant found in the Ras Responsive Element Binding Protein 1 gene (RREB1, gene ID: 6239) on chromosome 6 (Table 1). This m⁶A-SNP is part of an incomplete DRACH motif (TGA*AC). Near this m⁶A-SNP, a complete DRACH motif can also be observed at seven nucleotides downstream and two incomplete DRACH motifs localized at 14 and 25 nucleotides upstream in the RNA. Table 2 shows the occurrence for allele A, with a maximum frequency of 0.159 in EAS populations and the minimum frequency of 0.006 in African (AFR) populations. Moreover, this m⁶A-SNP alters a VDR motif, which could result in altered regulation of gene expression (Table 1).

The m⁶A-SNP rs76563149 (G/T) showed a cis-eQTL effect in different cells and tissues: fibroblast, adipose tissue, artery, oesophagus, lung, skeletal muscle, tibial nerve, stomach, testis, thyroid, and whole blood. This m⁶A-SNP appears in the 5’-UTR region of the Zinc Finger Protein 354A gene (ZNF354A, gene ID: 6940) on chromosome 5 (Table 1). Since this gene is oriented in the minus strand, the variants of this m⁶A-SNP in the mRNA will be C/A. Variant A of this m⁶A-SNP belongs to an incomplete DRACH motif (GAA*GC). In this case, the frequency of allele A is in the range of 0.081 to 0.368 for AFR and South Asian (SAS) populations, respectively (Table 2). Moreover, this m⁶A-SNP could be involved in gene expression regulation by altering a regulatory motif and the binding of 11 proteins (Table 1).

Finally, the m⁶A-SNP rs11614913 (C/T) showed cis-eQTL effect in different cells and tissues: fibroblast, tibial artery and nerve, skin, and whole blood. The rs11614913 is present in the microRNA 196a-2 gene (MIR196A2, gene ID: 406,973) between Hox10 and Hox9 genes on chromosome 12 (Table 1). This m⁶A-SNP is located near two complete DRACH motifs (at 4 and 19 nucleotides upstream in the RNA) and one incomplete DRACH motif (at 11 nucleotides upstream in the RNA), and, therefore, it could influence the m⁶A modification (Figure 2). In this case, the frequency of allele C is in the range of 0.458–0.860 for EAS and AFR populations, respectively (Table 2). The change from the U to the C variant does not alter the secondary structure of this miRNA (Figure 2). However, this m⁶A-SNP could also be involved in gene expression regulation by altering a regulatory motif.

Figure 2. A graphical representation of the miRNA structure obtained from the UNAFold Web Server (http://www.unafold.org/mfold/applications/rna-folding-form.php) [52]. The m⁶A-SNP rs11614913 (C/U) is found inside the red circle and it may influence the m⁶A modification of the complete DRACH motif (green boxes) and incomplete one (blue box). (a) Representation of the miRNA structure with the variant C for the m⁶A-SNP rs11614913. (b) Representation of the miRNA structure with the variant U for the m⁶A-SNP rs11614913.

Differential expression analysis

For the four genes, in which the m⁶A-SNPs showed cis-eQTL effect, we analysed the gene expression level in breast tissue. In the cases of the m⁶A-SNPs rs1801270, rs9379084, and rs76563149 that are present in the CDKN1A, RREB1, and ZNF354A genes, respectively, the gene expression levels were obtained from two gene expression microarray data sets from GEO database. In the case of the m⁶A-SNP rs11614913 present in the MIR196A2 gene, the gene expression level was obtained from the ENCORI database. We found that only ZNF354A and MIR196A2 present a differential gene expression. With respect to ZNF354A, differential expression was obtained in the GSE21422 study (p = 0.019) that analysed five healthy tissue samples and 14 carcinoma samples (9 for ductal carcinomas in situ and 5 for invasive ductal carcinomas). For MIR196A2, a significant difference in gene expression was found when comparing 1,085 samples of breast invasive carcinoma tissue with 104 normal tissues (p = 4.0e-14).

Although CDKN1A did not present significant differences in gene expression in the two datasets analysed, the study of Wei et al. (2015) [53] suggested otherwise. They analysed the gene and protein expression of CDKN1A in 65 breast tumour samples and their corresponding noncancerous adjacent tissues and their results showed that the differences in gene and protein expression of CDKN1A were statistically significant (p < 0.01). For RREB1, no statistical significance was found after analysing the data from GEO database, and we did not find any related study of gene expression.

Validation of the m⁶A-SNPs in patients with breast cancer

We determined the genotype of 35 BC patients for the three m⁶A-SNPs (rs1801270, rs76563149, and rs11614913) in order to calculate the genotype and allele frequencies. We compared these frequencies with data available in the 1000 Genomes database for females in the IBS population. In the group of patients with BC, the percentage of individuals with the A allele was 11.40% (AA = 5.70% and AC = 5.70%) for the m⁶A-SNP rs1801270 and 65.70% (AA = 14.30% and AC = 51.40%) for the m⁶A-SNP rs76563149 (Table 3). In the case of the m⁶A-SNP rs11614913, the percentage of individuals that presented the C allele, that is the variant associated with BC in different studies [54–58], was 80.00% (CC = 51.40% and CT = 28.60%). On the other hand, in the IBS population, 17.00% (AA = 0% and AC = 17.00%) and 45.30% (AA = 3.80% and AC = 41.50%) showed the A allele for the m⁶A-SNPs rs1801270 y rs76563149, respectively, while the 83.00% (CC = 43.40% and CT = 39.60%) had the C allele for the m⁶A-SNP rs11614913 (Table 3). No statistically significant differences were detected between the patients with BC and the IBS population when we compared the frequencies of the genotypes for the three m⁶A-SNPs (Table 3). However, when we compared allele frequencies, statistically significant differences were observed in the m⁶A-SNP rs76563149 (p = 0.029).

Table 3. Alleles and genotypes frequencies and χ² test for the patients with BC and female from the IBS population. We took values of p < 0.05 as significant (indicated in bold in the table below).

		Frequencies
m^6A-SNP/Gene	Alleles/Genotypes	Controls (n = 53)	BC patients (n = 35)	χ^2	P-value
rs1801270 CDKN1A	AA + AC	0.170	0.114	0.516	0.472
	CC	0.830	0.886
	A	0.085	0.086	<0.001	0.985
	C	0.915	0.914
rs76563149^a ZNF354A	TT + GT	0.453	0.657	3.536	0.060
	GG	0.547	0.343
	T	0.245	0.400	4.745	0.029
	G	0.755	0.600
rs11614913 MIR196A2	CC + CT	0.830	0.800	0.129	0.719
	TT	0.170	0.200
	C	0.642	0.657	0.045	0.832
	T	0.358	0.343

^aThe variants for this m⁶A-SNP in the mRNA are A or C.

Since the m⁶A-SNP rs76563149 located in the ZNF354A gene was the only one that presented significant differences, we proceeded to evaluate its expression in the tumour breast tissue and its corresponding adjacent normal tissue for our BC patients. A significant difference was observed between the tumour and normal breast tissue (p < 0.001), with a reduced expression in the tumour samples (0.785 ± 0.613 in tumour breast tissue and 1.095 ± 0.564 in normal breast tissue). Given that BC tumorigenesis/prognosis is affected significantly by menopausal status, we evaluated the ZNF354A expression considering this variable. No significant difference was found when compared the menopausal status of the patients (p = 0.522). Subsequently, we evaluated the ZNF354A expression taking into account the molecular subtypes. In this case, only the luminal B HER2+ subtype presented a significant difference (p = 0.001), with a reduced expression in the tumour breast tissue (Table 4 and Figure 3a-e). No significant differences were found when comparing the tumour subtypes with each other (p = 0.117) (Table 4 and Figure 3f).

Table 4. Average of relative gene expression for the m⁶A-SNP rs76563149 localized in ZNF354A gene with its standard deviation (SD) in the tumour breast tissue and its corresponding adjacent normal breast tissue and the p-values obtained for the Mann–Whitney U and Kruskal–Wallis H tests. We took p-values of <0.05 as significant, except in the Kruskal–Wallis H test, where a p-value of <0.005 was considered as significant (indicated in bold in the table below).

		Luminal A	Luminal B HER2-	Luminal B HER2+	HER2+	TNBC
Average ± SD	Normal tissue	1.016 ± 0.191	1.151 ± 0.645	1.032 ± 0.270	1.274 ± 1.161	1.037 ± 0.304
	Tumour tissue	0.924 ± 0.657	0.752 ± 0.167	0.492 ± 0.163	0.622 ± 0.468	1.156 ± 1.064
Mann–Whitney U test (p)		0.180	0.180	0.001	0.109	0.222
Kruskal–Wallis H test (p)		0.117

Figure 3. Relative expression levels of the ZNF354A gene. The fragment analysed contains the m⁶A-SNP rs76563149. Gene expression between each tumour subtype and its corresponding adjacent normal tissues: (a) TNBC. (b) Luminal A. (c) Luminal B HER2 +. (d) Luminal B HER2-. (e) HER2 +. (f) Gene expression between the different tumour subtypes. * P-value <0.05. The circles are outlier.

Discussion

Recently, m⁶A-SNPs have begun to be considered an important class of functional genetic variants that may provide new insights into the onset, development, and progression of diseases. In this regard, m⁶A-SNPs may be involved through multiple effects. Among others, they can affect the levels of m⁶A methylation, affect the binding of transcription factors and proteins, modify the secondary structure of mRNA, and affect the production of mature mRNA with correct biological function when they are found at splice sites [13,14]. For this reason, several studies have been carried out in order to identify m⁶A-SNPs and evaluate the role that they play in different diseases [13, 14, 16, 17, 19–27].

In the present study, we found 981 m⁶A-SNPs after comparing the list of SNPs related to BC with the list of m⁶A-SNPs. Of these, only four m⁶A-SNPs presented a cis-eQTL effect: rs1801270 in CDKN1A, rs76563149 in ZNF354A, rs9379084 in REEB1, and rs11614913 in MIR196A2, and all of them, except rs9379084, showed an altered gene expression in breast tissue when using data from GEO and ENCORI databases and Wei et al. (2015) [53]. These m⁶A-SNPs are different from those found by Xuan et al. (2021) [27], so our study allows us to expand the list of m⁶A-SNPs that could be involved in BC.

m⁶A-SNP rs9379084 present in the RREB1 gene did not show differences in gene expression after analysing the data from the GEO database. However, the protein encoded by this gene acts in different biological processes, such as DNA damage repair, cell growth and proliferation, cell differentiation, fat development, glucose balance fasting, zinc transport, and transcriptional regulation [59], which may play a role in cancer. In fact, RREB1 has been involved in the development of several cancers and other diseases [59], indicating that it may also be related to the development of BC although the effects of the methylation of this m⁶A-SNP would be involved in processes subsequent to gene transcription.

To understand the potential role of these three m⁶A-SNPs identified in the present study and their corresponding genes in BC, we investigated their biological functions. The m⁶A-SNP rs11614913 is in the MIR196A2 gene. MIR196A2 encodes two different mature miRNAs, miR-196a-5P and miR-196a-3P, from the same stem-loop and rs11614913 is in miR-196a-3P. This mature miRNA regulates the expression of multiple target genes, mainly by pairing to the 3′-UTR region of mRNA, and affecting either stability, translational inhibition, or degradation of mRNA. Indeed, under different conditions, this miRNA may function as an oncomiR [60]. We genotyped this m⁶A-SNP in our BC patients and we found no significant difference in relation to the IBS population (Table 3). These results are in agreement with those obtained by other authors in Caucasian populations [54, 61, 62]. However, analysis of Asian populations shows significant associations where variant C increases the BC risk [54–58, 60]. Therefore, we believe that it is important to evaluate the role of this m⁶A-SNP in different ethnic groups. On the other hand, as shown in Figure 2, the presence of the C or U variant does not alter the secondary structure of this miRNA, and since this m⁶A-SNP is near four m⁶A positions, future work should address how this m⁶A-SNP may affect the m⁶A modification and its influence on BC.

The m⁶A-SNP rs1801270 is present in the CDKN1A gene. The expression of CDKN1A is controlled by the tumour suppressor protein p53 and thus it links DNA damage to cell cycle arrest. The protein encoded by this gene is involved in multiple functions. On the one hand, the CDKN1A protein binds and inhibits the activity of the cyclin-dependent kinases, functioning as a regulator of cell cycle progression at G1 when the DNA is damaged [63,64]. On the other hand, this protein interacts with proliferating cell nuclear antigen (PCNA) and plays a regulatory role in S phase DNA replication and DNA damage repair. Moreover, it also plays a role in apoptosis [65]. Given the involvement of CDKN1A in these cellular processes, different works have studied its role in the development, progression, and therapy of cancer [66,67]. In this sense, in comparison with their corresponding noncancerous adjacent tissues, BC samples show an increased mRNA and protein expression of CDKN1A [53]. On the other hand, several studies have associated the m⁶A-SNP rs1801270 with different types of cancer and reported that the presence of the A allele is a risk factor [68–73]. However, the studies carried out on BC in the European population [74,75] did not confirm its association with variant A, whose frequency in this population is the lowest (Table 2). This corresponds to the experimental results obtained in the present study for patients with BC. In our patients, the presence of the A allele did not present significant differences with the IBS population. Therefore, this m⁶A-SNP should be evaluated in other populations where its frequency is higher, e.g. in the EAS, AFR, or AMR populations, in order to determine its involvement in BC. Despite not finding differences, we hypothesize that since m⁶A-SNP rs1801270 is a missense variant found in the coding region, this m⁶A-SNP could not only have an amino acid change but could also affect the RNA processing by altering the stability, nuclear export, transcription elongation, or splicing of the CDKN1A mRNA in BC [11, 13–15]. Moreover, this m⁶A-SNP could regulate the translation by reducing the kinetics of codon-anticodon pairing [15] or resolving secondary structures of mRNA [15,76].

The m⁶A-SNP rs76563149, located in the ZNF354A gene, belongs to the Krüppel‐associated box domain zinc finger proteins (KRAB-ZNF), one of the most abundant families of epigenetic repressors [77]. KRAB-ZFPs silence transcription by the deposition of repressive chromatin marks [77]. Moreover, KRAB-ZFPs may also facilitate DNA methylation adjacent to their binding motifs, especially in cells with stem cell phenotype [77]. To date, the role of rs76563149 remains unexplored. This m⁶A-SNP is part of a DRACH motif and, consequently, its presence could influence m⁶A methylation. We genotyped this m⁶A-SNP in our BC patient samples and our results showed a higher frequency of allele A in the patients than in the females of IBS population (Table 3). As a result, we carried out a gene expression analysis and found low levels of expression in tumour breast tissues when compared with their corresponding adjacent normal breast tissues. When we carried out the gene expression analysis taking into account the molecular subtypes, the luminal B HER2+ subtype was the only one that presented significant differences. These results suggest that m⁶A-SNP rs76563149 could alter the m⁶A methylation and it may affect the expression of ZNF354A. In this way, it could affect the splicing, nuclear export, and stability of mRNA, as well as the binding of reader proteins [76]. Indeed, it has been observed that this m⁶A-SNP could alter the binding site of transcription factors and affect the binding of 11 proteins (Table 1). Moreover, a decrease of ZNF354A expression in tumour BC tissues could influence cell proliferation, migration, invasion, and apoptosis, like other best-studied members of KRAB-ZNF family [78]. Further studies are needed to evaluate the role that ZNF354A may play in luminal B HER2+ subtype because no data was found in the literature. However, since this molecular subtype is characterized as oestrogen receptor-positive (ER+), we hypothesize that ZNF354A could act as a transcriptional regulatory factor that modulates the ER pathway, like ZNF496 [79,80]. This member of KRAB-ZNF family interacts with the DBD domain of ER type α (ER-α), inhibiting the binding of ER-α to the oestrogen response element, which selectively reduces the expression of ER-α target genes and inhibits the growth of ER-α-positive BC cells [79,80]. However, in the case of ZN354A the low expression would produce an opposite effect, increasing the growth of ER-α-positive BC cells. In addition, the effect of ZNF354A could be enhanced when the tumour is not only ER+ but also human epidermal growth factor receptor 2-positive (HER2+), as is the case of luminal B HER2 +.

While in this work we have focused on the m⁶A-SNPs that can influence gene expression and are statistically significant, the remaining m⁶A-SNPs that did not present these characteristics should not be ruled out. The m⁶A methylation could have more subtle effects that imply a fine modulation of gene expression and, in this case, the influence of methylation could not be observed by the analysis of eQTL. For this reason, more experimental studies should be carried out in order to confirm whether these 978 m⁶A-SNPs are involved in BC.

In conclusion, our in silico analysis identified three m⁶A-SNPs that may be associated with BC. These m⁶A-SNPs, rs11614913, rs1801270, and rs76563149 are in the MIR196A2, CDKN1A, and ZNF354A genes, respectively. After in vitro analysis, only the m⁶A-SNP rs76563149 shows significant differences. Future investigations of these m⁶A-SNPs should expand the study in different ethnic groups and increase the sample sizes to test their association with BC and elucidate their molecular function.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data Availability Statement

The data that support the findings of this study are available from the corresponding authors, MMP and FO, upon request.

Supplementary material

Supplemental data for this article can be accessed online at https://doi.org/10.1080/15592294.2022.2111137

Additional information

Funding

This work was supported by Basque Government (Grupos Consolidados with Grant No. IT1633-22 and predoctoral fellowship PRE_2019_2_0005 to Tamara Kleinbielen).