Candidate l‐methionine target piRNA regulatory networks analysis response to cocaine‐conditioned place preference in mice

Abstract Background Methionine has been proven to inhibit addictive behaviors of cocaine dependence. However, the mechanism of methionine response to cocaine CPP is unknown. Recent evidence highlights piRNAs to regulate genes via a miRNA‐like mechanism. Here, next‐generation sequencing is used to study mechanism on methionine response to drug‐induced behaviors though piRNA. Methods l‐methionine treatment cocaine CPP animal model was used to do non‐coding RNA sequencing. There were four groups to sequence: saline+saline (SS), MET+saline (MS), MET+cocaine (MC), and cocaine+saline. Combining mRNA sequencing data, the network and regulation of piRNA were analyzed with their corresponding mRNA and miRNA. Results Analysis of the piRNAome reveals that piRNAs inversely regulated their target mRNA genes. KEGG analysis of DE‐piRNA target mRNA genes were enriched in Morphine addiction, GABAergic synapse and Cholinergic synapse pathway. Furthermore, four significantly differential expressed genes Cacna2d3, Epha6, Nedd4l, and Vav2 were identified and regulated by piRNAs in the process of l‐methionine inhibits cocaine CPP. Thereinto, Vav2 was regulated by multiple DE piRNAs by sharing the common sequence: GTCTCTCCAGCCACCTT. Meanwhile, it was found that piRNA positively regulates miRNA and three genes Bcl3, Il20ra, and Insrr were identified and regulated by piRNA through miRNA. Conclusion The results showed that piRNA negatively regulated target mRNA genes and positively regulated target miRNA genes. Genes located in substance dependence, signal transduction and also nervous functions pathways were identified. When taken together, these data may explain the roles of l‐methionine in counteracting the effects of cocaine CPP via piRNAs.


INTRODUCTION
Addictions are a diverse set of common, complex diseases that are to some extent tied together by shared genetic and environmental etiological factors. Cocaine affects many neurochemical systems. Previous studies have further demonstrated that PFC-dependent processes, such as executive function, explicit learning, and memory, are damaged in animal models of cocaine addiction and in human cocaine abusers (Chao et al., 2014). Knowledge of genetic factors in etiology and treatment response to cocaine addiction may enable the individualization of prevention and treatment, as well as the identification of new therapeutic targets.
Non-coding RNAs (ncRNAs) have been gaining recognition for their involvement in genetic and epigenetic regulation (Cech & Steitz, 2014). A large body of studies have reported the distribution of piR-NAs in brains, including the hippocampus. Increasing evidence suggests that piRNAs play a critical role in both epigenetic and posttranscriptional silencing of transposons (Siomi et al., 2011;Thomson & Lin, 2009). piRNAs are shown to direct DNA methylation on a non-transposon locus to regulate genomic imprinting and influence synaptic plasticity of neuronal cells through DNA methylation on nontransposon genomic regions (Rajasethupathy et al., 2012;Watanabe et al., 2011). There is evidence that the Piwi/piRNA complex in the Aplysia brain was sensitive to serotonin modulation, regulates CREB2 and genes to control nervous system function via a miRNA-like mechanism that operates by imperfect base-pairing rules (Aravin et al., 2001;Rouget et al., 2010;L. Zuo et al., 2016). Hence, conducting piR-NAome analysis is very important for uncovering mechanism of drug abuse.
Epigenetic changes, specifically alterations in the pattern of DNA methylation, which can produce long lasting alterations in gene expression affecting behavior, have been verified by mounting reports (Maze & Nestler, 2011;Tsankova et al., 2007). L-methionine (MET) as a dietary methyl donor is the major driver of DNA methylation (N. Zhang, 2018).
MET has partitioned functions among protein synthesis, redox balance, polyamine biosynthesis, and the de novo pathway (also referred to as the methylation cycle or recycling pathway) where it is converted to S-adenosyl methionine (SAM), a principal methyl donor (Bolander-Gouaille & Bottiglieri, 2007;Sanderson et al., 2019). MET is involved in the synthesis of various neurotransmitters in the brain. It is considered to be the principal sulfur-containing amino acid in proteins, and play critical roles in nurture-based metabolism with brain development and functioning (Ringman & Coppola, 2013;Vuaden et al., 2012). Numerous evidence indicates that dietary methionine restriction is associated with increased longevity and decreased incidence of age-related disorders in mice and rats (Malloy et al., 2006;Miller et al., 2005;Zimmerman et al., 2003). More importantly, our initial findings and that of others suggest that methionine administration can inhibit addictive behaviors in rodent models of cocaine dependence (Tian et al., 2012(Tian et al., , 2016Wright et al., 2015). Currently, the specific mechanisms of methionine effects on drug-induced behaviors remain unclear.
Based on this evidence, the goal of this study was to identify the mechanisms of MET inhibitory effects on cocaine-induced cellular and behavioral changes through piRNAs. In the current study, it is the first time we identified Vav2 genes regulated by multiple piRNAs in the PFC upon cocaine treatment and reversed by L-methionine. The novel targets will be useful for further therapeutic intervention.

Animals housed, drug treatment and CPP test
Animals housed, drug treatment and CPP were as described before (Wright et al., 2015). Briefly, adult male C57/BL6 mice (20-30 g) were housed under a 12 h light/dark cycle. Mice were injected subcutaneously with 1 g/kg (6.6 mmol/kg) L-methionine twice a day for 10 consecutive days using the CPP procedure. During training, mice were injected with methionine 1 h before each behavioral experiment.
Mice were paired up for 8 days with the saline group receiving saline in both sides of the chambers, and drug groups were injected with cocaine (20 mg/kg, i.p., Qinghai pharmaceutical group co. LTD, china) and saline on one side and drug only on the opposite side (Romieu et al., 2008). After each manipulation, the mice were confined to the corresponding conditioning chambers for 30 min before being returned to their cages. On test day, place preference score (CPP score) was

RNA extraction
Two hours after the (day 10) CPP test, the animals were killed (being anesthetized with dry ice) and their mPFCs surgically excised. The tissue was stored in liquid nitrogen immediately after extraction and then transferred to the −80 • C freezer. Total RNA was extracted from the frozen tissues using the TIANamp DNA/RNA isolation kit (TIANGEN), including additional treatment with RNase-free DNase I (Ambion) for 30 minutes at 37 • C to remove contaminating DNA.

piRNA sequencing data
In this study, piRNA sequencing data was separated from ncRNA sequencing data. 1 µg of total RNA from a pool of 2 animals (500 ng each) was used to separate and recover 18-30 nt RNA segments by PAGE gel. Then the standard step was followed for generating sequencing library from the addition of poly-A. Sequencing was performed with single-end reads of 50 bp through the BGISEQ-500 platform at Beijing Genomics Institute (BGI; Shenzhen, China). Three biological replicates were sequenced for each group.

piRNA bioinformatics
Raw sequencing reads were filtered to remove low quality, adapter contaminated, shorter than 16 nt reads with FASTQ (Cock et al., 2010).
Clean reads were mapped to the reference genome by using AASRA (Tang et al., 2017). The results matching to miRNA, rRNAs and tRNAs were excluded. The remaining reads were aligned against piRBase (J. Wang et al., 2019) using Bowtie allowing one mismatch. In order to investigate the expression profiles of piRNAs, the frequency of piRNA counts were normalized to TPM (tags per million) using the following formula: normalized expression = actual miRNA read count /total clean read count × 10 6 .
Differentially expressed piRNAs between the paired groups were analyzed by using DEGseq (L. Wang et al., 2010). The p-values calculated for each gene were adjusted to Q-values for multiple testing corrections by two alternative strategies (Benjamini & Hochberg, 1995;Storey & Tibshirani, 2003). To improve the overall accuracy of DEGs results, a gene was defined as a DE-piRNA (differentially piRNA) when Q-values≤ .05.

Real-time PCR validation
Validation of DEGs was performed by quantitative reversetranscription PCR (RT-qPCR) using the Maxima SYBR Green qPCR Master Mix kit (Fermentas), according to the manufacturer's instructions, in an ABI Prism 7500 Sequence Detection System machine (Applied Biosystems Inc.). All real-time RT-qPCR data were normalized to SS expression using forward and reverse primers listed in Supporting Information Table S1.

Statistical analysis
Data are expressed as the mean ± SD. Statistical analysis of behavior data and qPCR was performed using an unpaired t test with two tailed distributions. The results were considered statistically significant when p < 0 .05.

Ethics approval
This study has been approved by the ethics committee of the Institute of Psychology, CAS.

Effects of methionine on cocaine-induced behavior
As shown in Figure 1, there were significant differences in cocaineinduced CPP among groups of saline and cocaine CPP with or without methionine treatment (p < 0.001). Meanwhile, methionine administration significantly attenuated cocaine-CPP behavior (CS vs. MC p < 0.01).

Characterization of ncRNA sequencing
To understand the molecular mechanisms underlying the effects of cocaine CPP and delineate the functional consequences of methionine treatment, we performed ncRNA sequencing of the extracted PFCs of mice from each of the four treatment groups: SS, CS, MS, and MC.
We obtained on an average 10 million clean reads for each sample.
There was no significant difference of the reads among each group

The response of MET on cocaine CPP
Behavior test above showed that MET attenuated behaviorally cocaine-CPP which induced the question whether MET clearly reversed cocaine-induced piRNA genes. We performed ncRNA sequence and first plotted the log 2 fold_change of piRNAs from CS, MS, and MC conditions in a heatmap. As expected, a significant number of genes exhibited fold changes in opposite directions in the MC mice compared with CS ( Figure 3a,b, p-value: 2.44e-14). This indicated that L-methionine reversed piRNAs induced by cocaine. When we compared piRNAs of CS and MS groups, there is a positive relationship between them (Figure 3a,c). Considering that L-methionine did not induce any abnormal behavior, piRNAs induced by L-methionine treatment cannot be considered to relate to cocaine response.

Broad effect of cocaine and L-methionine on piRNAome
To investigate the effect of piRNA of saline and MET response to drug, differentially expressed piRNAs analysis of the ncRNA sequencing data was performed. There were 143 differentially expressed piRNAs  Table S2).
In order to study the effects of combined cocaine and MET treatment on piRNAs, we applied the same analyses to find piRNAs of which the expression was altered by repeated cocaine treatment, with the condition that differential expression is inhibited by MET treatment.
We performed a differential expression analysis between MC and CS groups to obtain a list of DE-piRNAs affected by MET treatment in the context of cocaine treatment. A total of 139 DE-piRNAs were found upregulated and 119 DE-piRNAs down-regulated in the MC-treated mice when compared with that of CS-treated mice (Figure 3f).
Then, we examined how behavioral differences between the effects of cocaine, MET and cocaine combining MET might be explained partially by the differences of piRNAs. The Venn diagram in Figure 3g reveals that there is significant overlapping of DE-piRNAs between CS_SS and MC_CS (p < 1.468e-07; p < 8.852e-31) and between MS_SS and MC_CS (p < 1.987e-05; p < 5.435e-54) but not between CS_SS and MS_SS, which reinforce the speculation that L-methionine treatment was not related to cocaine response in Figure 2c

piRNAs target mRNA gene prediction
Recent evidence suggests that piRNAs, in a mechanism similar to miR-   (Figure 4d). There is a highly significant (p < 6.18e-09) overlap of piRNAs up-regulated by cocaine and down-regulated by MET with target genes in opposite direction, and also highly significant (p < 6.34e-06) overlap piRNAs down-regulated by cocaine and up-regulated by MET with target genes in opposite direction. Then we did pathway analysis with these target genes and plotted in Sankey diagram (Figure 4e). It showed that except Acyp2 gene, Cacna2d3, Epha6, Nedd4l and Vav2 were downregulated by MET via piRNA when induced by cocaine. And as for pathways these genes located such as mmu04024: cAMP signaling pathway, mmu04010: MAPK signaling pathway, mmu04360: Axon guidance and so on were reported related to drug dependence (Hope et al., 2007;Y. Wang et al., 2018).

MET acts on piRNAs to down-regulate genes located in addiction-related pathways induced by cocaine
It was noted that there were multiple piRNAs induced by MET down-regulated Vav2. Vav2 is an important gene located in multiple pathways report associated with drug addiction such as mmu04024: cAMP signaling pathway, mmu04062: Chemokine signaling pathway (Cui et al., 2014), mmu04510: Focal adhesion (Verma et al., 2019) and so on as a hub. Identifying recurring sequence patterns is very important in order to understand how piRNA binds and controls its target mRNA. We used clustalX (http://www.clustal.org/) to do sequence alignment of the piRNAs and showed the common sequence: GTCTCTCCAGCCACCTT as the target sequence works on (Figure 4f).

piRNAs target miRNA prediction
It is reported that miRNA can regulate the expression of protein-coding mRNA transcripts by binding to the 3′ untranslated region (3′ UTR) of target transcripts and blocking their translation (Kenny, 2014). There are few papers which studied a combined miRNA-piRNA regulation or network to detect diseases (Jain et al., 2019). In current paper, we tried to draw the regulatory network of how piRNA regulates protein-coding genes via miRNAs.
According to base pair complementarity rules, we first performed network by significantly expressed piRNAs and their potential target miRNAs from our ncRNA sequencing data with cytoscape (v3.8.2) (Figure 5a). Mostly, unlike an mRNA gene can be targeted by multiple piR-NAs, a miRNA can be targeted by only one or several piRNAs. This is probably because there are thousands and tons of piRNAs and miRNAs discovered in mice. Until now, the number of unique piRNA sequences in the mouse is over 68 million (Kim, 2019). Then we plotted the cor- Interestingly, all three genes were reported to be associated with drug addiction (Freeman et al., 2010;Heller et al., 2015;Przybycien-Szymanska et al., 2014). Especially, Bcl3 belongs to Bcl3-NFKB2 complex which was reported to be involved in mediating complex behaviors including learning and memory, stress responses and drug reward (Nennig & Schank, 2017). This indicates that the MET can reverse F I G U R E 6 Expression of target genes regulated by piRNAs which was acted by cocaine and MET as determined by qPCR. (a) Normalized mRNA expression quantified by RT-qPCR for genes whose expression is altered by cocaine treatment and reversed by MET: Vav2 and Nedd4l.

Validation of genes reversed by L-methionine
Two gene of Vav2 (Vav Guanine Nucleotide Exchange Factor 2) and Nedd4l (Neural Precursor Cell Expressed, Developmentally Down-Regulated 4-Like) were identified targeted and regulated by DE-piRNAs in the process of MET response of cocaine CPP. RT-qPCR assay confirmed that expression of these two tested genes was significantly altered by cocaine and reversed by MET, which coincided with the results obtained from RNA-sequencing ( Figure 6).

DISCUSSION
Until recently, knowledge of the impact of abused drugs on gene and protein expression in the brain was still limited. As a class of important small ncRNA, piRNA gained a growing concern. More and more studies focused on their correlation with diseases. In this study, we assessed piRNA expression and tried to unclose the mechanism of piRNA in the process cocaine exposed and MET reversing cocaine-CPP with ncRNA sequencing data. Of all sequencing data, 7.32% reads are piRNA. 86% of piRNAs located in repetitive elements: LINE1 (35%), LTR/ERVK (27%), and LTR/ERVL-MaLR (24%). It was noted that piRNA usually matches repetitive elements, and gene expression regulated via repetitive elements is possible.
It was reported that the Piwi/piRNA complex regulates some protein-coding genes in nervous system (Lee et al., 2011), whereas the mechanism of piRNA in the process of drug addiction is unknown. Our results showed that there was an inverse correlation between piRNA expression and its corresponding mRNA targets as was reported (Hashim et al., 2014). Pathway analysis showed that Morphine addic- So far, there are few papers that discuss the regulation between miRNA and piRNA. In this paper, according to base pair complementarity rule, for the first time we found piRNA targeted miRNA with a possible positive correlation. At the same time, three down-regulated DEGs were identified and regulated by piRNA through miRNA in the process of L-methionine response of cocaine CPP: Bcl3, Il20ra and Insrr, which were reported to have an expression significantly regulated when drug acted on. Thus, L-methionine treatment has an inverse response to cocaine administration via up-regulating piRNAs and miRNA is possible.
In summary, our results revealed novel molecular mechanisms of MET that inhibited the rewarding actions of cocaine in brain reward circuitries, and provided a theoretic support for the development of antiaddiction therapeutics based on MET.

CONCLUSIONS
Our results showed that piRNA negatively regulated target mRNA genes and positively regulated target miRNA genes. Genes located in substance dependence, signal transduction and nervous functions pathways were down-regulated: Cacna2d3, Epha6, Nedd4l and Vav2.
Thereinto, Vav2 was targeted by multiple piRNAs by sharing the common target sequence: GTCTCTCCAGCCACCTT. Taken together, these data may explain the roles of L-methionine in counteracting the effects of cocaine CPP via piRNAs.

ACKNOWLEDGMENTS
This work was supported by the National Natural Science Foundation of China (No. 61401459, 91132728, 31741062) and Key Laboratory of Mental Health, Institute of Psychology, Chinese Academy of Sciences.

CONFLICTS OF INTEREST
The authors have no conflicts of interest to declare.

DATA ACCESSIBILITY
Non-coding RNA-seq gene expression data are deposited in NCBI GEO: Series GSE146631.

DATA AVAILABILITY STATEMENT
Data that support the findings of this study are available from the corresponding author upon reasonable request.

AUTHORS CONTRIBUTION
Kunlin Zhang analyzed the data; Guanyu Ji analyzed the data and edited the article; Mei Zhao designed the study and interpreted the data; Yan Wang performed experiments, analyzed the data, and drafted the article.