Acupuncture Treatment Reverses Retinal Gene Expression Induced by Optic Nerve Injury via RNA Sequencing Analysis

Glaucoma and traumatic optic nerve crush (ONC) injury result in progressive loss of retinal ganglion cells (RGCs) and defects in visual function. In clinical trials of Traditional Chinese Medicine, acupuncture has been widely used for the treatment of ocular diseases. However, the molecular mechanisms of acupuncture treatment are still unclear. In this study, we used technique of RNA sequencing (RNA-seq) to study the effects of acupuncture treatment on retinal transcriptome after axotomy injury. RNA-seq results revealed that 436 genes including 31 transcription factors (TFs) were changed after injury, among them were many well-known neural degeneration related TFs such as Jun, Ddit3, Atf3, and Atf4. Interestingly, acupuncture treatment at acupoint GB20 (Fengchi) significantly reversed a series of differential expressed genes (DEGs) induced by optic nerve injury. While treatments at BL1 (Jingming) or GB20 sham control acupoint-GV16 (Fengfu), led to limited DEG reversal. In contrast, treatments at these two sites further enhanced the trend of DEG expression induced by axotomy injury. At last, retina immunostaining results revealed that only GB20 acupoint treatment increased RGC survival, in consistent with RNA-seq results. Therefore, our study first reported that acupuncture treatment regulated retinal transcriptome and reversed the gene expression induced by axotomy injury, and GB20 acupoint treatment increased RGC survival, which will provide novel therapeutic targets for treatment of ocular diseases.


INTRODUCTION
Chronic eye diseases such as glaucoma and optic neuritis result in progressive loss of retinal ganglion cells (RGC) and blindness finally (Shindler et al., 2008;Almasieh et al., 2012). Traumatic injury such as optic nerve crush (ONC) also cause serious RGC death (Sanchez-Migallon et al., 2016). Previous studies show that both chronic and traumatic optic nerve injury induce endoplasmic reticulum (ER) stress in RGC soma, while manipulation of ER stress pathway promoted the survival of RGC Huang et al., 2017). The induction of upregulation of pro-apoptotic TFs such as JUN and DDIT3 after injury leads to progressive RGC death (Huang et al., 2017;Syc-Mazurek et al., 2017a,b;Welsbie et al., 2017). Acupuncture has been widely used to treat many kinds of pains such as neck and shoulder pain, lower back pain worldwide (MacPherson et al., 2017;Koh et al., 2018). Acupuncture has also been used for treatment of ocular diseases, including glaucoma, age-related macular degeneration (AMD), retinitis pigmentosa, etc (Jiao, 2011;Xu et al., 2012;Law and Li, 2013;Xu and Peng, 2015). In Traditional Chinese Medicine, GB20 and BL1 are the common acupoints selected for acupuncture treatment of ocular diseases (Qin et al., 2015;Xu and Peng, 2015). Although acupuncture is widely used for ocular therapy, the effect of acupuncture treatment on regulation of retinal transcriptome especially those changed by ONC injury is totally unknown. RNA sequencing (RNA-seq) is a new technique effective in identifying numerous genes regulated by specific treatment. In this study, we will use ONC crush injury mouse model and RNA-seq technique to analyze retinal samples with and without optic nerve injury and screen for the upregulated and downregulated TFs 2 days after ONC injury, and also identified those injury-induced genes reversed after acupuncture treatment at two different acupoints.

Expression Analysis
To investigate expression differentiation, we extract RNA from intact retina (WT), axotomized retina (ONC), and axotomized retina with acupuncture treatment at acupoint GB20 (ONC-F) or BL1 (ONC-J), 2 days after ONC. The actual and schematic acupuncture sites were shown in Figures 1A,B. RNA integrity of each sample was assessed by Bioanalyzer 2100 system. Sample total reads ranged from 40.8 to 69.5 million, and mapping rate ranged from 95.4 to 96.6%.
We then used feature Counts v1.5.0-p3 to calculate expected number of Fragments Per Kilobase of transcript sequence per Millions base pairs sequenced (FPKM) of each sample, which represent relative gene expression abundance. The boxplot result showed that the overall distribution of the FPKM values were consistent among samples, suggesting that the RNA-seq data were reproducible (Figure 2A). Principal component analysis (PCA) showed that wild-type (WT) control formed clear cluster apart from injury groups, and ONC with or without acupuncture treatment formed one big cluster, ONC with sham treatment was separated from acupuncture treatment groups ( Figure 2B). Differentially expressed genes (DEG) cluster analysis showed that acupuncture treatment groups were clustered together and separated from ONC control, and all the ONC groups with or without acupuncture treatments were separated from WT control group ( Figure 2C).

Differentially Expressed Genes After ONC
To identify candidate genes regulated by axon injury 2 days after axotomy, we performed differential expression analysis using the DESeq2 R package (1.16.1). Genes with | log 2 FoldChange| > 0 and adjusted p-value < 0.05 were assigned as differentially expressed. We then identified total 436 differentially expressed genes (DEG), with 191 genes upregulated and 245 genes downregulated ( Figure 3A). Enriched gene ontology analysis results showed that the most upregulated gene categories were cell adhesions and junctions ( Figure 3B); while the most downregulated gene categories were: axon, postsynapse and transmembrane transports ( Figure 3C). The 20 most significantly upregulated and downregulated DEGs after ONC are shown in Table 1.

Reversal of Injury Induced DEGs by Acupuncture Treatments
To study the effect of acupuncture treatment after ONC injury, we performed acupuncture treatments 3 times after ONC at acupoint GB20 (Fengchi), BL1 (Jingming) or sham point GV16 (Fengfu) (Figure 1B), respectively. RNA-seq analysis showed that among the injury induced up-regulated DEGs, 20 DEGs were reversed by GB20 treatment with only 6 DEGs further up-regulated ( Figure 4A), the reversal DEG number was 11 in BL1 and 14 in GV16 treatments, however, both BL1 and GV16 treatments further up-regulated even more DEGs, which was 24 in BL1 and 52 in GV16 (Figures 4B,C). Among the injury induced down-regulated DEGs, 21 DEGs were reversed by GB20 treatment without further down-regulated DEGs (Figure 4D), the reversal DEG number was 19 in BL1 with only 2 DEGs further down-regulated ( Figure 4E), however, in GV16 treatments, the reversal number is 32 vs. 42 further down-regulated DEGs ( Figure 4F). The average FPKM values of the 436 DEGs of all sample groups and the reversal percentages of all the DEGs after BL1, GB20, and GV16 treatments were listed in Supplementary Table S1.

Validation DEGs With qRT-PCR
To verify the reliability of RNA-seq result, we selected four important TFs reported to be upregulated after axon injury (Jun, Ddit3, Atf3, and Sox11) (Hu et al., 2012;Fagoe et al., 2015;Holland et al., 2016;Huang et al., 2017;Lu et al., 2017;Norsworthy et al., 2017;Struebing et al., 2017;Li et al., 2018;Wong et al., 2018) to perform qRT-PCR test. Results showed that all the four TFs were significantly upregulated after ONC injury (Figures 5A-D), consistent with previous reports. We then    also test the DEGs expression reversed by acupuncture treatment. Results showed that expression of Penk (Preproenkephalin) and Mt3 (Metallothionein-3) were significantly decreased after ONC injury, Penk was reversed by BL1 (ONC-J) treatment and Mt3 was reversed by GB20 (ONC-F) treatment, while sham treatment GV16 (ONC-S) had no effect on expression reversal of both candidate DEGs (Figures 5E,F). Expression of Klf6 (Kruppellike factor 6) and Nav1 (neuron navigator 1) were significantly increased after ONC injury and only GB20 treatment reversed and decreased the expression of candidate genes, while BL1 treatment had no effect on expression change and GV16 sham treatment further increased the expression of candidate genes (Figures 5G,H). All these qRT-PCR results were consistent with RNA-seq analysis and demonstrated that acupuncture treatment regulated retinal transcriptome after ONC injury.

Acupuncture Treatment Promote RGC Survival
Finally, to test how acupuncture treatment affected actual survival of RCGs after ONC, we collected and immunostained the retinas after ONC alone and ONC with acupuncture treatment at acupoints BL1, GB20, and sham control acupoint GV16. Retinas were co-stained with RGC specific marker-RBPMS and neuron specific marker-Tuj1. Results showed that GB20 acupoint treatment slightly increased RGC survival after optic nerve injury. 27% of RCGs survived in GB20 acupuncture condition, while only 18% survived under ONC alone. BL1 and GV16 had no significant effect on RGC survival, which was 18 and 22%, respectively (Figure 6). Thus, ONC with acupoint GB20 treatment resulted in 9% more RCGs surviving compared to ONC alone.

DISCUSSION
In this study, we performed RNA-seq analysis of retinal transcriptome 2 days after ONC. Axotomy injury induced both upregulation and downregulation of more than eight hundred genes, including many pro-apoptotic molecules and ER stress pathway molecules. Our previous studies showed that both traumatic and chronic Frontiers in Integrative Neuroscience | www.frontiersin.org neuron injury resulted in activation of PERK-eIF2α-ATF4-CHOP pathway or downregulating activity of PI3K-AKT-mTOR pathways, while manipulated these pathways and reversed the injury induced changes resulted in neuroprotection and increased RGC survival (Park et al., 2008;Miao et al., 2016;Yang et al., 2016;Huang et al., 2017Huang et al., , 2019. DLK-JNK-JUN pathway is another important pathway activated by axon injury and knockdown this pathway also resulted in neuroprotection (Kim et al., 2016;Syc-Mazurek et al., 2017b).
In Traditional Chinese Medicine, acupuncture has been broadly used for treatment of ocular diseases. GB20 and BL1 are the most common acupoints for ocular disease treatment. However, BL1 acupoint is too close to the eyeball and it's difficult to handle during clinical practice, and it is the same in mouse acupuncture treatment (Figure 1); GB20 acupoint has many muscles surrounded and it is easy to perform needling at this acupoint. RNA-seq analysis revealed that acupuncture treatment at GB20 reversed the expression of 41 genes and only 6 genes were regulated toward the same direction induced by ONC injury (Figures 4A,D). In BL1 acupoint treatment, 43 genes were reversed, and 13 genes were regulated toward the same direction induced by ONC injury (Figures 4B,E). Although sham treatment at GV16 also reversed 46 genes, 94 genes were regulated toward the same direction induced by ONC injury (Figures 4C,F). As in qRT-PCR results, sham treatment further increased Klf6 and Nav1 expression which were already increased by ONC injury (Figures 5G,H). Therefore, only GB20 treatment had the best overall reversal effect on DEGs expression, and GB20's sham control GV16 treatment had no overall reversal effect at all, in contrast, sham treatment further enhanced the expression states induced by ONC injury. What is more, the promotion of RGC survival by GB20 treatment confirmed the neuroprotection effect of DEG reversal (Figure 6B).
Among these reversal DEG candidates, Mt3 plays an important role in ocular neovascularization and its deficiency will exacerbate retinal degeneration (Tsuruma et al., 2012; FIGURE 4 | Reversal of DEGs by acupuncture treatment after ONC injury. (A-C) Volcano map shows the reversal of up-regulated DEGs induced by ONC injury after acupuncture treatments at acupoint BL1, GB20, and GV16 (Sham control site). (D-F) Volcano map shows the reversal of down-regulated DEGs induced by ONC injury after acupuncture treatments at acupoint BL1, GB20, and GV16. p value < 0.05. Choi et al., 2013); Fbn1 (fibrillin 1) expression is required for eye development and its mutation is associated with macular degeneration (Hubmacher et al., 2014;Ratnapriya et al., 2014), suggesting that acupuncture treatment may affect vascularization processes in the retina, and promote neuroprotective outcomes after ONC injury. Intriguingly, the site of acupoint GB20 is at the back of the neck, far from the site of eyes, while needling and stimulation at this point surprisingly has effect on regulation of retina gene expression, which explains why many clinical practices choose GB20 site to cure ocular diseases. At last, it is possible that acupuncture treatments relay electrical signal via dorsal root ganglion (DRG) and spinal cord, finally regulate ocular blood flow since many groups reported that acupuncture treatment improved eye blood flow in openangle glaucoma patients (Leszczynska et al., 2018;Vanzini and Gallamini, 2018). It is reported that acupuncture treatment also regulated expression of nerve growth factor and brain-derived neurotrophic factor in retina (Pagani et al., 2006) which may regulate retinal gene expression and provide neuroprotection.
Since GB20 treatment slightly increase RGC survival, next we will investigate the detail mechanism under GB20's involvement in RGC protection.

Animals
We perform experiments in 6 weeks old C57BL/6 mice. All animal procedures were performed in accordance with the National Institute of Health guidelines. The protocol was approved by the Animal Care and Use Committee of Peking University.

Optic Nerve Crush
Mice were anesthetized by xylazine and ketamine based on their body weight (0.01 mg xylazine/g + 0.08 mg ketamine/g). Optic nerves of both sides were sequentially exposed intraorbitally and crushed with a jeweler's forceps (Dumont #5; Fine Science Tools) for 2 s approximately 0.5 mm behind the eyeball. Care was taken not to damage the underlying ophthalmic artery. Erythromycin Eye Ointment was applied to protect the cornea after surgery.

Acupuncture Treatment
After optic nerves crush injury, mice were anesthetized by xylazine and ketamine based on their body weight (0.01 mg xylazine/g + 0.08 mg ketamine/g) before acupuncture treatment at acupoint GB20 or BL1 at both sides, respectively. The depth of the acupuncture needling is around 2mm. The duration of acupuncture treatment was 20 min. We then gave another two acupuncture treatments every 24 h.

RNA Preparation
Mice were randomly divided into four groups (3 mice/group). Experiments were repeated for 3 times. Briefly, in each replicate, mice were sacrificed 2 days after ONC injury, and retinas were dissected out in HBSS buffer (Cellgro) immediately. Retinas were then homogenized with TRIzol Reagent (Thermo Fisher Scientific), and total RNA was extracted from the homogenized mixture according to the reagent instructions. RNA purity was checked using the NanoPhotometer R spectrophotometer (IMPLEN). RNA concentration was measured using Qubit R RNA Assay Kit in Qubit R 2.0 Flurometer (Life Technologies). RNA integrity was assessed using the RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies).

Library Preparation and Sequencing
About 3000 ng total RNA generated from each group was used for RNA-seq, which was done at Novogene, Inc. Briefly, RNA samples from three biological replicates went through mRNA purification with poly-T oligo-attached magnetic beads. Sequencing libraries were generated using NEBNext R UltraTM RNA Library Prep Kit for Illumina R (NEB) following manufacturer's recommendations and index codes were added to attribute sequences to each sample. The library fragments were purified with AMPure XP system (Beckman Coulter) for cDNA fragments of preferentially 250∼300 bp in length. Library quality was assessed on the Agilent Bioanalyzer 2100 system. The clustering of the index-coded samples was performed on a cBot Cluster Generation System using TruSeq PE Cluster Kit v3-cBot-HS (Illumia) according to the manufacturer's instructions. After cluster generation, the library preparations were sequenced on an Illumina Hiseq platform and 125 bp/150 bp paired-end reads were generated.

Gene Expression Analysis
The RNA-seq reads were aligned to the reference genome using Hisat2 v2.0.5. FeatureCounts v1.5.0-p3 was used to count the reads numbers mapped to each gene. And then FPKM of each gene was calculated based on the length of the gene and reads count mapped to this gene, which normalizes gene expression by considering the effect of sequencing depth and gene transcript length at the same time. Differential expression analysis was performed using the DESeq2 R package (1.16.1). The resulting p-value were adjusted using the Benjamini and Hochberg's approach for controlling the false discovery rate (less than 0.05). Genes with an adjusted p-value < 0.05 (detected by DESeq2) were considered to be differentially expressed.

Quantitative Real-Time PCR (qRT-PCR)
Twelve RNA samples (3 samples each group) were examined with qRT-PCR. 1 µg total RNA from each sample was reversetranscribed into cDNA using SuperScript III (Invitrogen). All qRT-PCR assays were performed in duplicate. Reactions were carried out with SYBR R Premix Ex Taq TM kit (TaKaRa) and performed on LightCycler R 480 System (Roche). For comparison of relative gene expression, we analyzed qRT-PCR data by Ct method and normalized value to endogenous GAPDH control.

Immunostaining of Flat-Mount Retina
Retinas were dissected out from 4% PFA fixed eyes and washed extensively in PBS before blocking in staining buffer (10% normal goat serum and 2% Triton X-100 in PBS) for 30 min. Mouse neuronal class ß-III tubulin (clone Tuj1, 1:200; Biolegend), rabbit RBPMS (GTX118619, 1:400, GeneTex) were diluted in the same staining buffer. Floating retinas were incubated with primary antibodies overnight at 4 • C and washed three times for 30 min each with PBS. Secondary antibodies (Cy2 and Cy3conjugated) were then applied (1:200; Jackson ImmunoResearch) and incubated for 1 h at room temperature. Retinas were again washed three times for 30 min each with PBS before a cover slip was attached with Fluoromount-G (Southernbiotech).

Counting Surviving RGCs
For RGC counting, whole-mount retinas were immunostained with RBPMS antibody, and 4-6 fields were randomly sampled from peripheral regions of each retina. The percentage of RGC survival was calculated as the ratio of surviving RGC numbers in injured eyes compared to contralateral uninjured eyes.

Statistical Analysis
Data are represented as means ± SEM. qRT-PCR data was analyzed with Student's t-test, and P-value < 0.05 was considered as statistically significant. The raw data and GEO accession number for this study are as follows: GSE131486, link: https://www.ncbi.nlm.nih.gov/geo/query/acc. cgi?acc=GSE131486.

DATA AVAILABILITY STATEMENT
The datasets generated for this study can be accessed on NCBI website via GEO accession number: GSE131486.

ETHICS STATEMENT
All animal procedures were performed in accordance with the National Institutes of Health guidelines.