GPX8 deficiency–induced oxidative stress reprogrammed m6A epitranscriptome of oral cancer cells

ABSTRACT

Glutathione peroxidase 8 (GPX8) is a key regulator of redox homoeostasis. Whether its antioxidant activity participates in the regulation of m⁶A modification is a crucial issue, which has important application value in cancer treatment. In this study, MeRIP-seq was used to explore the characteristics of transcriptome-wide m⁶A modification in GPX8-deficient oral cancer cells. Oxidative stress caused by the lack of GPX8 resulted in 1,279 hyper- and 2,287 hypo-methylated m⁶A peaks and 2,036 differentially expressed genes in GPX8-KO cells. Twenty-eight differentially expressed genes were related to the cell response to oxidative stress, and half of them changed their m⁶A modification. In GPX8-KO cells, m⁶A regulators IGF2BP2 and IGF2BP3 were upregulated, while FTO, RBM15, VIRMA, ZC3H13, and YTHDC2 were downregulated. After H₂O₂ treatment, the expression changes of RBM15, IGF2BP2, and IGF2BP3 were further enhanced. These data indicated that GPX8-mediated redox homoeostasis regulated m⁶A modification, thereby affecting the expression and function of downstream genes. This study highlights the possible significance of GPX8 and the corresponding m⁶A regulatory or regulated genes as novel targets for antioxidant intervention in cancer therapy.

Key policy highlights

• Lack of GPX8 caused oxidative stress of oral cancer cells.

• Oxidative stress induced by GPX8 deficiency reprogrammed m⁶A epitranscriptome.

• GPX8 deletion–caused oxidative stress regulated expression of m⁶A regulatory genes.

• m⁶A modification of antioxidant genes is the adaptive response of cells to oxidative stress.

KEYWORDS:

• glutathione peroxidase 8 (GPX8)
• oxidative stress
• m⁶A modification
• reactive oxygen species (ROS)
• m⁶A regulatory genes

Introduction

Glutathione peroxidase (GPX) family proteins play an important role in maintaining cell redox balance and normal cell function [1]. GPX8 is the most recently identified member of this family. Its peroxidase activity and endoplasmic reticulum (ER) localization can prevent the leakage of H₂O₂ from ER and maintain the ER redox control [2]. GPX8 has a broad range of biological functions. It was found to be a cellular substrate of the hepatitis C virus NS3-4A protease [3]. GPX8 can regulate Ca²⁺ storage and flux in ER and its expression was related to the concentration of Ca²⁺ in the ER and the flux of Ca²⁺ in cytoplasm and mitochondria [4]. GPX8 protected insulin-secreting INS-1E beta-cells against lipotoxicity by improving the ER antioxidative capacity [5]. It was reported that GPX8 is a negative regulator of caspase 4/11, which can prevent against colitis [6]. GPX8 deficiency impacted on the lipid composition of cancer cell microsomal membranes [7]. GPX8 was also identified as a key gene involved in the spermatogenesis in patients with cryptorchidism [8].

There is increasing evidence that GPX8 is dysregulated in various tumours. Existing studies have shown that GPX8 is related to tumour progression and prognosis. Bioinformatics analysis revealed that the expression of GPX8 was positively correlated with progression and poor prognosis in patients with gastric cancer [9,10], pancreatic cancer [11], colon cancer, and lung cancer [12,13]. GPX8 can also inhibit the apoptosis of tumour cells and promote their migration and invasion by regulating epithelial properties [13]. Elevated GPX8 activated Wnt signalling pathway to promote the proliferation, migration, and invasion of gastric cancer cells [14]. Histone deacetylase inhibitors inhibited the expression of GPX8, which makes hepatocellular carcinoma sensitive to ER stress and apoptosis through oxidative stress [15].

Oxidative stress is characterized by the excessive production of reactive oxygen species (ROS), which is widely recognized as a key factor in many pathophysiological processes and cancer development [16]. ROS has dual biological properties, playing an anti-tumour or tumour promoting role in different tumours [16,17]. Elevated levels of ROS are commonly observed in cancer cells and generally play a tumour promoting role, but too high level of ROS is toxic to cancer cells [18]. Normal cells adopt several mechanisms to maintain intracellular ROS levels and overall redox homoeostasis. They can express antioxidant enzymes such as superoxide dismutase, catalase, and glutathione peroxidase and produce nonenzymatic antioxidants such as glutathione and thioredoxin, to protect them from ROS damage [19]. In order to adapt to relatively high levels of ROS, cancer cells must enhance their antioxidant capacity to neutralize the cytotoxicity of excessive ROS [17]. Antioxidant protection treatment converts ROS into less reactive species, which can neutralize the harmful effect of ROS [16].

Conventional chemotherapy and radiotherapy are closely related to oxidative stress and can cause ROS-mediated DNA damage and apoptosis. For example, the generation of ROS contributes to the cytotoxicity of 5-fluorouracil, and antioxidant treatment facilitates the drug resistance of tumour cells to 5-fluorouracil [20,21]. The arsenic-mediated elevation of ROS induces apoptosis and pyroptosis in a variety of cancer cells [22,23]. Ionizing radiation-mediated elevation of ROS induces oxidative damage and ferroptosis [24–26]. Excessive levels of ROS induced by ionizing radiation in radiotherapy will destroy the redox homoeostasis, leading to oxidative stress that may result in cell death. On the contrary, tumour cells remove excess ROS by activating endogenous antioxidant enzymes, thereby generating radioresistance [27]. Although antioxidant therapy is a potential strategy for ROS-induced cancer, antioxidant treatment may increase the risk of some cancers and promote their progression [16,28]. Therefore, the regulation of oxidative stress is related to the efficacy of tumour therapy.

N⁶-methyladenosine (m⁶A) is the most common eukaryotic mRNA modification, which widely regulates RNA transcription, maturation, translation, and metabolism, thus affecting various physiological and pathological processes, including oxidative stress and tumorigenesis. The modification of m⁶A is a dynamic and reversible process, which is coordinated by methyltransferase/m⁶A writer and demethylases/m⁶A eraser. The m⁶A writers include methyltransferase like-3 (METTL3)/METTL14/METTL16, Wilms tumour 1-associated protein (WTAP), RNA binding motif protein 15 (RBM15), RBM15B, Vir-Like m⁶A methyltransferase associated (VIRMA), and zinc finger CCCH-type containing 13 (ZC3H13), which transfer methyl to the target sites. The m⁶A erasers mainly include fat mass and obesity-associated protein (FTO), AlkB homolog 5 (ALKBH5) and ALKBH3, which can remove m⁶A from the modified sites of RNAs. The m⁶A reader proteins recognize m⁶A modification of RNAs, mediating downstream biological functions. Several classes of proteins act as m⁶A readers, including YT521-B homology (YTH) domain-containing proteins YTHDC1/2 and YTHDF1/2/3, insulin-like growth factor 2 mRNA binding protein 1 (IGF2BP1)/IGF2BP2/IGF2BP3, heterogeneous nuclear ribonucleoproteins HNRNPA2B1 and HNRNPC, and eukaryotic initiation factor 3 (eIF3).

The m⁶A modification is closely related to oxidative stress. The cross-talk between oxidative stress and m⁶A modification and its significance in the occurrence and development of tumours are very complicated. First, ROS can dynamically regulate the expression and activity of m⁶A regulators, thereby changing cellular m⁶A level [29,30]. On the other hand, m⁶A modification of the genes related oxidative stress may regulate their own expression, thus affecting the balance of oxidation and antioxidation, as well as consequent occurrence and progression of tumours [31,32]. Therefore, understanding the interaction mechanism between oxidative stress and m⁶A methylation in tumours is of great significance for tumour therapy. In a previous study, we found that GPX8 was upregulated in oral cancer cells. However, it is not clear whether its antioxidant activity involves the regulation of m⁶A modification in cells. In this study, we profiled the transcriptome-wide m⁶A methylome in GPX8-deficient oral cancer cells through MeRIP-seq and measured the expression of major m⁶A regulatory genes in order to understand the effect of oxidative stress induced by antioxidant gene deletion on m⁶A modification. This will help to evaluate the cross-talk between oxidative stress and m⁶A modification and find novel targets for tumour therapy.

Materials and methods

Construct and generation of GPX8 knockout cell line

Oral squamous cell carcinoma cell SCC-9 was provided by Shanghai Guandao Biological Engineering Co., Ltd. (Sgdbio, Shanghai, China) [33,34]. GPX8-KO SCC-9 cells was derived from SCC-9 cells, in which GPX8 gene was edited by clustered regularly interspaced short palindromic repeats (CRISPR)/CRISPR associated protein 9 (Cas9) system. Briefly, according to the sequence GPX8 gene, two sgRNAs targeting 5’ end of GPX8 open reading frame (aactgcaaatacctttgctc and tttagcgggtaagctgcaag) were designed and cloned into lentiCRISPR-v2 vector. The 293T cells were incubated with constructs and lentiviral packaging plasmids psPAX2 and pMD2.G for 66 h. The resulting supernatant containing lentiviral particles was filtered with 0.45 µm filter and concentrated. 2.2 × 10⁵ SCC-9 cells in a well of 6-well plate were infected with 500 µL of concentrated virus. After 48 h of infection, the cells were screened in the medium with 0.3 μg/mL puromycin. The concentration of puromycin-resistant cells was adjusted to 5 cells/mL, and the cells were seeded into a 96-well plate at 100 µl per well to isolate a single clone. The knockout of GPX8 gene in single clone was verified by DNA sequencing and Western blot.

Cells and cell culture

SCC-9 and GPX8-KO SCC-9 cells were cultured in Dulbecco’s Modified Eagle’s Medium (DMEM) supplemented with 10% FBS (GIBCO, Australia), 100 U/mL penicillin G, and 100 μg/mL streptomycin. All cells were cultured under humidified conditions with 5% CO₂ at 37°C. About 5 × 10⁷ cells were collected for MeRIP-seq and RNA sequencing (RNA-seq).

MeRIP-seq and bioinformatics analysis

Transcriptome-wide methylated RNA immunoprecipitation sequencing (MeRIP-seq) and RNA sequencing (RNA-seq) were performed as previously described [35–37]. Briefly, the total RNA was extracted using TRIzol reagent (Invitrogen, Carlsbad, CA, USA). PolyA RNA was purified using Dynabeads Oligo (dT)25–61005 (Thermo Fisher, CA, USA) and cleaved into about 100nt fragments using Magnesium RNA Fragmentation Module (NEB, USA). One portion of the RNA fragments was used as input and the other portion was immunoprecipitated with m⁶A-specific antibody (Synaptic Systems, Germany) to enrich m⁶A-methylated RNA fragments. The RNA sequence library was prepared and purified, and the paired-end sequencing (PE150) of libraries was performed on Illumina Novaseq™ 6000 platform (LC-Bio Technology CO., Ltd., Hangzhou, China), following the manufacturer’s protocol.

Low quality reads and adaptors were removed from the raw data using fastp tool (https://github.com/OpenGene/fastp) [38]. The quality of IP and input sequences was verified using FastQC (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/) and RseQC (http://rseqc.sourceforge.net/) [39,40]. HISAT2 (http://daehwankimlab.github.io/hisat2) was used to map clean reads to the reference genome Homo sapiens (Version: v101) [41]. m⁶A peak calling and differentially methylated peaks were analysed by exomePeak2 of R package (https://bioconductor.org/packages/release/bioc/html/exomePeak2.html) [42], and R package ANNOVAR (http://www.openbioinformatics.org/annovar/) [43] was used to annotate the peaks. MEME (http://meme-suite.org) [44] and HOMER (http://homer.ucsd.edu/homer/motif) were used for motif finding. StringTie (https://ccb.jhu.edu/software/stringtie) [45] was used to analyse the expression level of transcripts and genes through calculating FPKM, and the differentially expressed transcripts and genes were identified using R package edgeR (https://bioconductor.org/packages/edgeR) [46]. The differentially methylated coding genes and differentially expressed genes were analysed by Gene Ontology (GO) and Kyoto Encyclopaedia of Genes and Genomes (KEGG), respectively [47].

Real-time reverse transcription-polymerase chain reaction (RT-PCR)

Total RNA was isolated using TRIzol reagent (Invitrogen, Carlsbad, CA, USA) according to the manufacturer’s instructions. RNA was treated with gDNA Eraser to remove genomic DNA and reverse-transcribed using PrimeScript RT Enzyme Mix I (PrimeScriptTM RT reagent Kit with gDNA Eraser, TaKaRa, Shiga, Japan) and random hexamer primers. The RT-PCR was carried out in 20 μL reaction volumes in triplicate with SYBR® Premix Ex Taq II (Tli RNaseH Plus, TaKaRa, Shiga, Japan) and the ABI PRISM®7900 system (ABI). The threshold cycles and relative expression levels were calculated with 2^−ΔΔCt. The primers used to detect 19 major m⁶A regulatory genes have been described in a previous study [37]. To investigate the effects of oxidative stress, cells were treated with 100 μM hydrogen peroxide for 12 h or 24 h before RNA isolation.

ROS assay

2‘,7’-dichlorofluorescein diacetate (DCFDA) cellular ROS detection assay kit (ab113851, Abcam, USA) was used to detect the ROS level in wild-type and GPX8-KO SCC-9 cells. The cells were seeded in 35 mm cell culture dishes and maintained in the DMEM medium for 24 h. Cells were washed twice with 1× buffer and incubated with 200 μl 25 μM DCFDA working solution at 37 ºC in the dark for 45 min. After two washings, the cells were examined under a fluorescence microscope. For flow cytometry analysis, cells were trypsinized and incubated with 25 μM DCFDA at 37°C for 30 min. After washing twice, the cells were suspended in PBS and analysed by CytoFlex (Beckman Courtier) at 488 nm.

Statistical analysis

According to the criteria |log₂FC| ≥ 1 and P < 0.05, the differentially methylated peaks of MeRIP-seq and differentially expressed genes of RNA-seq were determined using exomePeak2 [42] and edgeR [46] in R package, respectively. The GraphPad Prism 8 was used to analyse RT-PCR data, and two tailed unpaired Student's t-test was used to calculate P values.

Results

GPX8 deficiency increased the production of cellular ROS

In order to explore the role of GPX8 in maintaining cell redox homoeostasis, we designed two gdRNAs and edited GPX8 gene in SCC-9 oral cancer cells using CRISPR-Cas9 technology. Single colonies of GPX8-edited cells were screened and sequenced. In a monoclonal cell line, a guanine nucleotide was inserted at codon 5 of the coding sequence of GPX8, causing a frameshift mutation in the open reading frame (ORF) of GPX8 and premature termination of translation at codon 43 (Figure 1a, b). Western blot conformed that there was no GPX8 proteins in the GPX8-edited cell line, which was named GPX8-KO SCC-9 cells (Figure 1c). When the cells were stained with DCFDA, it was found that the ROS level in GPX8-deficient SCC-9 cells was significantly higher than that in wild-type SCC-9 cells (Figure 1d, e), indicating that the lack of GPX8 led to oxidative stress in cells.

Figure 1. The deletion of GPX8 induced oxidative stress in cells. (a) Editing of the GPX8 gene resulted in the insertion of a guanine nucleotide at codon 5 of the GPX8 coding sequence. (b) The edited GPX8 gene caused a frameshift mutation and premature termination of translation at codon 43. (c) GPX8 expression was not detected in GPX8-KO SCC-9 cells. (d, e) Compared with wild-type SCC-9 cells, GPX8-deficient SCC-9 cells increased ROS production. **P < 0.01.

Loss of GPX8 reprogrammed m⁶A epitranscriptome in oral cancer cells

To explore the effect of oxidative stress induced by GPX8 deletion on the m⁶A modification, we performed MeRIP-seq and RNA-seq of GPX8-deficient SCC-9 and SCC-9 cells (Table S1). After removing the reads with adapters, low-quality bases, and undetermined bases in the raw data, more than 90% of the clean reads from the IP and input samples can be mapped to gene exons in the reference genome. A total of 45,108 and 43,608 m⁶A peaks were identified in SCC-9 and GPX8-deficient SCC-9 cells, respectively. The two cell lines shared similar characteristics in the distribution and density of m⁶A peaks that were highly enriched in 3`UTR and stop codon regions (Figure 2a). Typical conserved m⁶A motifs were enriched in some m⁶A peak sequences (Figure 2b, c).

Figure 2. Distribution of m⁶A peaks across the mRNA transcripts and the representative m⁶A motifs enriched. (a) The m⁶A peaks were highly enriched in 3`UTR and stop codon regions. (b) The m⁶A motifs with typical conserved sequence in SCC-9 cells. (c) The m⁶A motifs with typical conserved sequence in GPX8-KO SCC-9 cells. NC, SCC-9 cells; KO, GPX8-KO SCC-9 cells.

Then, we compared the difference of m⁶A peaks between two cell lines and found 1,279 hyper-methylated and 2,287 hypo-methylated m⁶A peaks in GPX8-KO SCC-9 cells compared with SCC-9 cells (|log₂FC|≥1.0 and P < 0.05). The 20 genes with the most significant m⁶A modification change are listed in Table 1. The GO function and KEGG pathway enrichment of the differentially methylated mRNA was performed to explore the biological significance of m⁶A modification. These genes were enriched into GO terms such as protein binding and KEGG pathways such as ubiquitin mediated proteolysis (Figure 3a, b). Further analysis indicated that many genes listed under GO term of cellular response to oxidative stress had m⁶A modification change (Table 2).

Figure 3. GO function and KEGG pathway enrichment of differentially methylated mRNA. (a) Top 20 significantly enriched GO terms. (b) Top 20 significantly enriched KEGG pathways.

Table 1. The 20 coding genes with the most significantly altered m⁶A peaks in GPX8-KO SCC-9 cells compared with SCC-9 cells.

Hyper-methylated		Hypo-methylated
Genes	Peak region	Genes	Peak region
TCAF1, TRPM8 channel associated factor 1	exonic	CDADC1, cytidine and dCMP deaminase domain containing 1	5’ UTR
EIF3C, eukaryotic translation initiation factor 3 subunit C	exonic	CXCL8, C-X-C motif chemokine ligand 8	exonic, 3’ UTR
RGPD6, RANBP2 like and GRIP domain containing 6	exonic, 3’ UTR	HSF2BP, heat shock transcription factor 2 binding protein	3’ UTR
SERF1B, small EDRK-rich factor 1B	3’ UTR	SP8, Sp8 transcription factor	5’ UTR
INHBE, inhibin subunit beta E	3’ UTR	SLC31A1, solute carrier family 31 member 1	5’ UTR
KLHL3, kelch like family member 3	3’ UTR	PBLD, phenazine biosynthesis like protein domain containing	3’ UTR
ROBO3, roundabout guidance receptor 3	5’ UTR	TANGO6, transport and golgi organization 6 homolog	intronic
CLEC18B, C-type lectin domain family 18 member B	3’ UTR	TCEANC2, transcription elongation factor A N-terminal and central domain containing 2	3’ UTR
RGPD2, RANBP2 like and GRIP domain containing 2	3’ UTR	MORF4L1, mortality factor 4 like 1	5’ UTR
CTAGE15, CTAGE family member 15	exonic	PKP4, plakophilin 4	5’ UTR

Table 2. Differential expression of oxidative stress-related genes in GPX8-KO SCC-9 cells.

Genes	log2FC	m^6A change
NCF2, neutrophil cytosolic factor 2	7.14	-
PTGS2, prostaglandin-endoperoxide synthase 2	4.70	-
STC2, stanniocalcin 2	3.84	-
ETV5, ETS variant transcription factor 5	2.80	3’ UTR, exonic/down
SLC7A11, solute carrier family 7 member 11	2.58	3’ UTR/down
NOX5, NADPH oxidase 5	2.16	-
PRNP, prion protein	1.88	3’ UTR/up
CYBB, cytochrome b-245 beta chain	1.83	-
RBPMS, RNA binding protein, mRNA processing factor	1.77	3’ UTR/up
SOD2, superoxide dismutase 2	1.76	3‘UTR, intronic/down; 3’ UTR/up
MCTP1, multiple C2 and transmembrane domain containing 1	1.76	-
ZC3H12A, zinc finger CCCH-type containing 12A	1.70	-
BCL2, BCL2 apoptosis regulator	1.65	-
DHRS2, dehydrogenase/reductase 2	1.33	-
NCOA7, nuclear receptor coactivator 7	1.30	3’ UTR, exonic/down
ATF4, activating transcription factor 4	1.26	5’ UTR/down
SLC1A1, solute carrier family 1 member 1	1.15	-
RBM11, RNA binding motif protein 11	1.10	-
SESN2, sestrin 2	1.10	-
CYBA, cytochrome b-245 alpha chain	−1.79	3‘UTR/up; 5’ UTR/down
PPARGC1A, PPARG coactivator 1 alpha	−1.73	3’ UTR/down
NQO1, NAD(P)H quinone dehydrogenase 1	−1.46	3‘UTR/up; 5’ UTR/down
PRKD1, protein kinase D1	−1.40	3’ UTR/up
GPX8, glutathione peroxidase 8	−1.40	5’ UTR/down
HSPA1A, heat shock protein family A (Hsp70) member 1A	−1.15	-
IDH1, isocitrate dehydrogenase (NADP(+)) 1	−1.10	3‘UTR/up; 5’ UTR/down
NR4A2, nuclear receptor subfamily 4 group A member 2	−1.10	3’ UTR/down
ENDOG, endonuclease G	−1.09	-

Change of m⁶A modification in GPX8-KO SCC-9 cells may regulate gene expression

It seems that the loss of GPX8 has no significant effect on the transcriptional activity of the whole cell (Figure 4a). When RNA-seq data were analysed for differentially expressed genes, 1,123 genes were significantly upregulated and 913 genes were significantly downregulated genes in GPX8-KO SCC-9 cells (|log₂FC|≥1.0 and P < 0.05; Figure 4b). Among them, 28 genes are involved in the cellular response to oxidative stress (Table 2). The top 20 and top 100 genes differentially expressed between GPX8-KO SCC-9 cells and SCC-9 cells are shown in Table 3 and Fig S1.

Figure 4. Transcriptional activity of the whole cell and differentially expressed genes. (a) GPX8-KO SCC-9 and SCC-9 cells have similar overall transcriptional activity. (b) Volcano plots of differentially expressed genes between GPX8-KO SCC-9 and SCC-9 cells. |log₂FC|≥1.0 and P < 0.05. KO, GPX8-KO SCC-9 cells; NC, SCC-9 cells.

Table 3. The top 20 differentially expressed coding genes in GPX8-KO SCC-9 cells compared with SCC-9 cells.

Upregulated	log2FC	Downregulated	log2FC
CXCL10, C-X-C motif chemokine ligand 10	∞	RANBP3L, RAN binding protein 3 like	−7.49
CCL20, C-C motif chemokine ligand 20	∞	C2CD4A, C2 calcium dependent domain containing 4A	−5.10
IFNL1, interferon lambda 1	∞	VAV3, vav guanine nucleotide exchange factor 3	−5.02
CXCL1, C-X-C motif chemokine ligand 1	9.62	LGSN, lengsin, lens protein with glutamine synthetase domain	−4.48
IL13RA2, interleukin 13 receptor subunit alpha 2	8.28	NOG, noggin	−4.40
CCL5, C-C motif chemokine ligand 5	8.17	FRMD3, FERM domain containing 3	−4.00
OASL, 2’−5’-oligoadenylate synthetase like	8.14	PRG4, proteoglycan 4	−3.90
CXCL11, C-X-C motif chemokine ligand 11	7.51	SLC12A3, solute carrier family 12 member 3	−3.80
CCND1, cyclin D1	7.46	PLK2, polo like kinase 2	−3.70
IL33, interleukin 33	7.40	IGFBP5, insulin like growth factor binding protein 5	−3.63

Note: FPKM>2 in up-regulated genes.

We further jointly analysed the MeRIP-seq and RNA-seq data and identified 509 upregulated and 453 downregulated mRNAs in GPX8-KO SCC-9 cells. These genes were either m⁶A hyper-methylated or m⁶A hypo-methylated (Figure 5, Table 4). Even a single mRNA, such as serpin family E member 1 (SERPINE1) or DAB adaptor protein 2 (DAB2), can be modified by hypermethylation and hypomethylation at different sites at the same time (Table 4). The data indicated that the expression changes of many genes in GPX8-KO SCC-9 cells appeared to be related to the m⁶A modification, although the regulation of mRNA m⁶A methylation on its own expression needed further evaluation. The GO function and KEGG pathway enrichment of these genes were further analysed (Figure 6).

Figure 5. Distribution of differentially expressed genes with differential m⁶A peaks. Hyper-up, m⁶A peak upregulated and mRNA expression upregulated; Hyper-down, m⁶A peak upregulated and mRNA expression downregulated; Hypo-up, m⁶A peak downregulated and mRNA expression upregulated; Hypo-down, m⁶A peak downregulated and mRNA expression downregulated.

Figure 6. GO function and KEGG pathway enrichment of differentially expressed genes with differential m⁶A peaks. (A) Top 20 significantly enriched GO terms. (B) Top 20 significantly enriched KEGG pathways.

Table 4. The first 10 differentially expressed genes with differential m⁶A peaks in GPX8-KO SCC-9 cells compared with SCC-9 cells.

Hypo-methylated genes	log2FC	Hyper-methylated genes	log2FC
Upregulation of expression
CXCL1, C-X-C motif chemokine ligand 1	9.62	CCL5, C-C motif chemokine ligand 5	8.17
DNER, delta/notch like EGF repeat containing	7.05	SLC6A15, solute carrier family 6 member 15	5.92
IFIT2, interferon induced protein with tetratricopeptide repeats 2	6.79	ALCAM, activated leukocyte cell adhesion molecule	5.23
FN1, fibronectin 1	6.56	RAET1L, retinoic acid early transcript 1 L	5.03
RSAD2, radical S-adenosyl methionine domain containing 2	6.30	INHBA, inhibin subunit beta A	4.55
IL6, interleukin 6	5.93	ANTXR2, ANTXR cell adhesion molecule 2	4.54
IFIT1, interferon induced protein with tetratricopeptide repeats 1	5.92	NRP1, neuropilin 1	4.37
SP8, Sp8 transcription factor	5.64	SERPINE1, serpin family E member 1	4.29
CXCL8, C-X-C motif chemokine ligand 8	5.45	DUSP6, dual specificity phosphatase 6	3.83
ADM2, adrenomedullin 2	4.77	CFB, complement factor B	3.74
Downregulation of expression
LGSN, lengsin, lens protein with glutamine synthetase domain	−4.48	FRMD3, FERM domain containing 3	−4.00
PLK2, polo like kinase 2	−3.70	PRG4, proteoglycan 4	−3.90
IGFBP5, insulin like growth factor binding protein 5	−3.63	GLP2R, glucagon like peptide 2 receptor	−3.55
DAB2, DAB adaptor protein 2	−2.71	SFRP1, secreted frizzled related protein 1	−3.47
TSHZ2, teashirt zinc finger homeobox 2	−2.53	KCNH5, potassium voltage-gated channel subfamily H member 5	−3.14
MTARC1, mitochondrial amidoxime reducing component 1	−2.50	MVD, mevalonate diphosphate decarboxylase	−2.97
STEAP4, STEAP4 metalloreductase	−2.33	GATD3B, glutamine amidotransferase like class 1 domain containing 3B	−2.80
ATP1B1, ATPase Na+/K+ transporting subunit beta 1	−2.29	DAB2, DAB adaptor protein 2	−2.71
CDK10, cyclin dependent kinase 10	−2.19	H3C10, H3 clustered histone 10	−2.49
SPIRE2, spire type actin nucleation factor 2	−2.18	H3C8, H3 clustered histone 8	−2.47

Note: FPKM>2 in up-regulated genes.

GPX8 deficiency–induced oxidative stress dysregulated expression of m⁶A regulatory genes

Since RNA m⁶A modification is dynamically regulated by m⁶A regulatory genes, such as m⁶A writer and eraser, the expression of these genes was further investigated. By analysing our sequencing data, we found that insulin-like growth factor 2 mRNA binding protein 2 (IGF2BP2) (log₂FC, 1.32) and IGF2BP3 (log₂FC, 3.49) were upregulated and fat mass and obesity-associated protein (FTO) (log₂FC, −1.26) was downregulated in GPX8-KO SCC-9 cells compared with SCC-9 cells (p < 0.05). These were confirmed when the expression of 19 m⁶A regulatory genes was further measured by RT-PCR (P < 0.01, Figure 7a). In addition, the expression of RNA binding motif protein 15 (RBM15), Vir-Like m6A methyltransferase associated (VIRMA), zinc finger CCCH-type containing 13 (ZC3H13), and YTH domain-containing 2 (YTHDC2) decreased significantly in GPX8-KO SCC-9 cells (P < 0.01). To some extent, methyltransferase like-3 (METTL3), RNA binding motif protein 15B (RBM15B), heterogeneous nuclear ribonucleoprotein A2B1 (HNRNPA2B1) and heterogeneous nuclear ribonucleoprotein C (HNRNPC) were downregulated in GPX8 deficient cells (0.01< P < 0.05).

Figure 7. Expression of m⁶A regulatory genes detected by real-time RT-PCR. (A) Expression of m⁶A regulators in GPX8-KO SCC-9 and SCC-9 cells. (B) Expression of m⁶A regulators after hydrogen peroxide treatment within the indicated time. All bars indicate the mean ± SD of three independent experiments. *P < 0.05, **P < 0.01.

To study the effects of oxidative stress on the expression of m⁶A regulatory genes, the cells were treated with hydrogen peroxide for 12 h or 24 h, and several m⁶A regulatory genes were detected by RT-PCR. We found that the expression of RBM15, IGF2BP2, and IGF2BP3 further decreased or increased in GPX8-KO SCC-9 cells after 24 h of treatment (P < 0.01, Figure 7B). The expression of FTO and YTHDC2 was also downregulated to some extent (0.01< P < 0.05). Hydrogen peroxide treatment induced relatively smaller changes in the expression of these genes in SCC-9 cells.

Discussion

GPX8 is a recently discovered member of the GPX family, which is upregulated in many tumours. The significance of its antioxidant capacity in tumours and its cross-talk with epigenetic modification of m⁶A has not been well understood. In this study, we analysed the transcriptome-wide m⁶A modification through MeRIP-seq in GPX8-deficient oral cancer cells and investigated the expression of major m⁶A regulatory genes by RT-PCR. We found that the loss of GPX8 resulted in over production of ROS in cells. Thousands of genes in GPX8-KO SCC-9 cells, including nearly half of the differentially expressed genes, have changed their m⁶A modification status, and many of them are involved in the cell’s response to oxidative stress. GPX8 deficiency also regulated the expression of several m⁶A regulatory genes, and hydrogen peroxide stimulation led to further changes in the expression of these genes.

We can see the m⁶A peaks of both hyper-methylation and hypo-methylation in GPX8-deficient cells, and the m⁶A peak distribution of GPX8-deficient cells and wild-type cells is similar, which indicates that the lack of GPX8 will not affect the overall m⁶A methylation level and m⁶A modification characteristics. Changes in m⁶A modification may vary by gene. Hyper- or hypo-methylation can occur at multiple sites of a single transcript.

We found that in GPX8-KO SCC-9 cells, many genes changed their m⁶A modification, and some of them showed differential expression, which needs to be evaluated in the future. First, we need to confirm the change of m⁶A modification through various experimental methods. It is reported that some genes identified here are regulated by m⁶A methylation. For example, m⁶A methylation of eukaryotic translation initiation factor 3 subunit C (EIF3C) mRNA dependent on WTAP facilitates prostate cancer metastasis [48]. Secondly, we are not sure whether the difference of m⁶A modification will lead to the change of gene expression level. Among 2036 (1123 + 913) differentially expressed genes in this study, 962 (509 + 453) genes changed the state of m⁶A modification, which supports the regulatory effect of m⁶A modification on gene expression. The mechanism by which some key genes regulate their expression through m⁶A modification is worthy studying.

The deletion of GPX8 may regulate gene expression by affecting mRNA m⁶A modification. In turn, these genes may function to adapt to cellular oxidative stress, which is caused by the loss of GPX8. We must point out that many genes with differential m⁶A modification or differential expression in GPX8-KO SCC-9 cells do not seem to be directly related to oxidative stress. Their potential association with oxidative stress and the significance of these changes in GPX8 deficient cells remain to be clarified. However, we found that 28 genes (Table 2) related to the cellular response to oxidative stress were differentially expressed in GPX8-KO SCC-9 cells. Half of them changed their m⁶A levels. This indicates that GPX8 loss mediated oxidative stress can modulate the expression of genes, including those regulating redox homoeostasis.

The significance of cross-talk between m⁶A modification and oxidative stress is an important issue in many biological processes [19]. m⁶A modifications can modulate cellular ROS levels, and in turn, ROS signalling also plays a regulatory role in m⁶A modifications. ROS significantly induced global mRNA m⁶A levels by regulating ALKBH5 post-translational modification [49]. Hypoxia can induce ROS production [50]. ALKBH5 reduced the overall m⁶A level in response of hypoxia [51]. Hypoxia-induced breast cancer stem cell phenotype through ALKBH5-mediated m⁶A demethylation of NANOG mRNA [52]. YTHDF1 played an important role in hypoxic adaptation and pathogenesis of non-small cell lung cancer [53]. Hypoxia can induce SUMOylation of YTHDF2, which promotes mRNA degradation and cancer progression by increasing its binding affinity with m⁶A modified mRNAs [29]. Hypoxia blocked ferroptosis of hepatocellular carcinoma via suppression of METTL14 triggered YTHDF2-dependent silencing of SLC7A11 [54]. Some studies also described the role of m⁶A modification in regulating the redox balance of cells. For examples, several m⁶A-related lncRNAs are associated with oxidative stress in oral cancer and can predict the prognosis of oral squamous cell carcinoma [55]. The demethylase FTO promoted oxidative stress and apoptosis of ovarian cancer cells to inhibit tumour growth in nude mice [56]. It can also inhibit oxidative stress by mediating m⁶A demethylation of Nrf2 to reduce cerebral ischaemia/reperfusion injury [57]. YTHDC2 inhibited lung adenocarcinoma tumorigenesis by suppressing SLC7A11-dependent antioxidant function [58].

In fact, it has not been reported that the antioxidant activity of GPX8 is related to m⁶A modification. In this study, based on the following findings, we concluded that the oxidative stress caused by GPX8 deficiency changed m⁶A modification. First, GPX8 deficiency increased ROS level in cells. Next, the hydrogen peroxide stimulation in GPX8-KO SCC-9 cells did change the expression of several m⁶A regulatory genes (Figure 7B). This indirectly supports the change of m⁶A modification mediated by oxidative stress, but the regulatory effect of these genes on m⁶A modification is still unclear. Finally, the global change of m⁶A modification was observed in GPX8-KO cells. In the future, it will be interesting to explore the direct correlation between oxidative stress induced by GPX8 deficiency and m⁶A modification. It should be noted that the biological function and pathological significance of GPX8 are largely attributed to its antioxidant capacity. So far, there is no evidence that GPX8 acts as an m⁶A regulatory gene.

In conclusion, our results suggest that GPX8 lack-induced oxidative stress reprograms m⁶A epitranscriptome of oral cancer cells. It is supposed that oxidative stress modulates the expression of m⁶A regulatory genes, which in turn leads to transcriptome-wide changes of m⁶A modification, and consequentially affects the expression and function of downstream genes to adapt to oxidative stress. The study highlights the potential value of GPX8 and the related m⁶A regulatory or regulated genes as novel targets for antioxidant intervention in cancer treatment.

Author contributions

DY, WZ, and SC contributed to the conception and design of the study. XC, LY, LZ, LC, YH, CW, and JW collected samples and performed experiments. All authors participated in the writing of the manuscript. All authors have read and agreed to the published version of the manuscript.

Disclosure statement

No potential conflict of interest was reported by the authors.

Data availability statement

All data generated or analysed during this study are available from the corresponding author upon reasonable request (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?) with GEO accession number GSE224718.

Supplementary material

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

Additional information

Funding

This study was funded by the National Natural Science Foundation of China (No. 81873711, No. 82073378 and No. 31670788) and by the Open Fund of Guangdong Key Laboratory of Pharmaceutical Functional Genes (No. 2020B1212060031).