To investigate how glia respond to neurotransmitters, we developed a protocol analyzing gene expression in adult fly glia (Figure 1A). Adult fly brains expressing a membranous GFP (UAS-mCD8GFP) under the control of a pan-glial driver repo-GAL4 (repo > mCD8GFP) were dissected and incubated with different types of neurotransmitters (10 mM) for 30 min. The neurotransmitters we used include: Glu, DA, ATP, ACh, and GABA. GFP-positive adult fly glia were purified by the fluorescence activated cell sorting (FACS) and total RNAs were extracted from these glia for further sequencing analysis. Expression changes of genes of interest were verified again by the reverse transcription quantitative PCR (RT-qPCR). Due to the heterogeneity in the size and shape of glial cells, adult fly brain cells were suspended into single cells and sorted by the criteria whether the cells express GFP (Figure 1B). Interestingly, the sorting plot revealed that 8–10% of total cells were GFP-positive, correlating with the reported number of glia in adult fly brains (Croset et al., 2018; Brunet Avalos et al., 2019; Park et al., 2022). After removal of cell debris, single GFP-positive glial cells were visualized by confocal microscopy (Figure 1C).

By RNA-sequencing analysis, we were able to analyze the glial response at the transcription level to each type of neurotransmitters; the differences in individual gene expression were compared. A total of 720 million reads of sequencing samples were obtained, with each transcriptome containing about 0.4 million reads in three biological replicates. A total of 16,107 genes were detected including 2,631 uncharacterized genes (annotated in CG number). Total 46,048 expressed transcripts including 18,763 new transcripts were identified. RSEM software was used to quantitatively analyze the expression level of transcripts, and TPM conversion was performed to homogenize the sequencing length and depth of transcripts, so to obtain standardized expression level of transcripts. According to the Violin plot showing the sample expression distribution, the gene expression level in each transcriptome is similar, ensuring the same amount of information collected from different transcriptomes (Figure 1D). In addition, about 9,000 genes were identified in each transcriptome as revealed by the sample expression distribution Upset plot, with over 7,800 overlapping genes identified in all six groups and the rest expressing in different groups (Figure 1E). Sample distance heatmap (Figure 1F) and PCA analysis (Figure 1G) revealed a strong correlation between the bio-replicated samples for each group treated with different neurotransmitters. These results suggest that the difference between replicates of each group is relatively large, whereas among replicates themselves is small. The RNA-sequencing results are therefore reliable, consistent, and qualified.

Next, genes differentially expressed among these groups were analyzed. Gene set analysis was carried out and the Upset plot revealed the total number of differentially expressed genes (DEGs) in adult fly glia upon neurotransmitter addition. A total of 2,259 DEGs were identified in Glu-treated glia, whereas only 371 DEGs were detected in ATP-treated glia (Figure 2A). Among these identified DEGs, 1,171 upregulated and 1,088 downregulated genes were identified in Glu-treated glia, whereas only 80 upregulated and 291 downregulated genes were identified in ATP-treated glia (Figure 2B). As more upregulated genes were identified in the Glu and DA groups, more downregulated genes were detected in the ACh, GABA, and ATP groups (Figure 2B). These results suggest that neurotransmitters have differential effects on the glial gene expression. A heatmap of gene cluster analysis was drawn to compare the expression level of individual gene in each group (Figure 2C). Z score was used to standardize the gene expression: when the score is lower than the average value, z is negative; otherwise, it is positive. As shown in the heatmap, genes in the top panel (green label on the left) were mostly upregulated in the DA-treated glia compared to the control group, whereas the same set of genes were mostly downregulated in the GABA-, ACh-, or Glu-treated glia. Genes in the middle panel (blue label on the left) were mostly downregulated in glia treated with all five types of neurotransmitters. At last, genes in the bottom panel (red label on the left) exhibited a low expression in the control group, yet mostly upregulated in the Glu-treated glia (Figure 2C). Consistently, changes in some of these genes involved in signal transduction, metabolism, and function of glial cells were also verified by RT-qPCR in our hand (Supplementary Figures 1, 2).

On the other hand, DEGs categorized by different analyses revealed further information on the function of affected genes (Supplementary Figures 3–5). For instance, COG classification showed that most genes affected in all five groups of glia were poorly characterized (Group S), revealing unrecognized glial cell functions that remain to be investigated. GO annotation and KEGG analysis showed similar trends in the number of genes involved in different cellular processes or pathways for all five types of neurotransmitter-treated glia. Considering that neurotransmitters likely regulate glia in similar mechanisms, it is feasible that no dramatic change in one single type of cellular events is detected.

Despite that glia are not excitable cells, they utilize similar repertoires of synthesis enzymes, receptors, and transporters when transmitting the cellular signals in response to neurotransmitters (Fellin and Carmignoto, 2004; Fiacco and McCarthy, 2006; Hartl et al., 2007; de Tredern et al., 2021). Next, we analyzed the expression change in a list of genes characterized for their function in the synthesis, transport, and receptor-based signaling of neurotransmitters. Genes selected are the synthesis enzymes, transporters, and receptors for the five types of neurotransmitters and they are listed in Table 1.

As shown by the Volcano plot and heatmap in Figure 3, significant changes in the expression of genes listed in Table 1 were uncovered. We used “+” or “−” symbol for indicating genes with increased expression (the adjusted p-value < 0.05, log2FC > 1) or decreased expression (the adjusted p-value < 0.05, log2FC < −1), respectively in Table 1. A color tile along with the number designation was shown in the upper-right corner of the heatmap in Figure 3 for indicating the difference in gene expression. By default, a three color scale is used with blue indicating low expression values (2–4-fold lower), white indicating similarly expressed genes, and red representing highly expressed genes (2–4-fold higher). Genes with expression differentially affected by the neurotransmitters were selected and shown in the Volcano plot. Overall, neurotransmitters failed to evoke significant change in the expression of synthesis enzymes and triggered more change in the receptor expression. Among the five types of neurotransmitters, Glu activate the expression of Table 1 genes more significantly. Whereas Glu mainly evoked upregulation in receptor expression, GABA rarely caused any increase, correlating with its inhibitory function (Figures 3D, E).

Specifically, ACh increased the expression of glutamate decarboxylase 1 (GAD1) and nicotinic acetylcholine receptor alpha 1, 5-7 (nAChRα1, 5-7), yet decreased the expression of Glutamate receptor IIA (GluRIIA) (Control vs. ACh, Figure 3A). ATP also decreased the GluRIIA expression, and increased the expression of nAChRα1, 5-7 as potently as ACh (Control vs. ATP, Figure 3B). DA upregulated the GABA transporter GAT and the Excitatory amino acid transporter 1 (Eaat1), both are transporters for neurotransmitter trafficking (Control vs. DA, Figure 3C). Whereas GABA did not elicit a dramatic change in the overall gene expression, Glu elicited increased expression in nAChRα1, 5-7, Dopamine/Ecdysteroid receptor (DopEcR), metabotropic GABA-B receptor subtype 1 (GABA-B-R1), metabotropic Glutamate Receptor (mGluR), and Glutamate receptor IB (GluRIB) (Figure 3E).

It has been widely recognized that glia are part of the tripartite synapse and respond to neuronal signals by calcium wave propagation, a functional way to relay signals to adjacent or distal glia via gap junctions (Kielian and Esen, 2004; Deitmer and Rose, 2010; Verkhratsky, 2010; De Bock et al., 2013; Mu et al., 2019; Spray et al., 2019). To further investigate whether neurotransmitters trigger change in genes involved in these aspects, a list of genes encoding channels, transporters, and gap junctions were compiled and analyzed (Table 2). Volcano plot and heatmap were also illustrated in Figure 4, with changes in the expression level of these genes described similarly in Table 1 and Figure 3. Interestingly, expressions of most genes in these categories were upregulated upon ATP, DA, or Glu addition, whereas some inhibitory effects were seen in the group of ACh- and GABA-treated glia. These results suggest that in general neurotransmitters trigger glial activation.

Specifically, the expression of Sarco/endoplasmic reticulum Ca(2+)-ATPase (SERCA) was upregulated in the Glu-treated glia, suggesting that Glu triggers glial calcium activity. The expression of nervana 3 (nrv3) was upregulated in ATP- and Glu-treated glia, whereas nervana 2 (nrv2) was more highly expressed in the DA and GABA group. For gap junction genes, Innexin 2 (Inx2) and Innexin 3 (Inx3) were both upregulated in the DA-treated glia, yet downregulated in ACh- and GABA-treated glia. Moreover, genes encoding channels such as Glutamate-gated chloride channel subunit alpha (GluClalpha) and Ae2 in Glu group, and transient receptor potential-like (trpl), sodium potassium chloride cotransporter (NKCC) and SLC22A family member (SLC22A) in the DA group were more highly expressed. On the other hand, the expressions of Shawn and NKCC were downregulated in the DA- or ACh- and GABA-treated glia, respectively (Figure 4C). Taken together, these results suggest that neurotransmitters trigger changes in the expression of genes functioning in glial activity.

Accumulating evidence has indicated that metabolic changes in neurons associate with brain dysfunction (Mlody et al., 2016; Diaz-Garcia et al., 2017; Onkar et al., 2020; Rickman et al., 2020; Bonvento and Bolanos, 2021; Li et al., 2021). On the other hand, not a lot has been explored for how glia react to neuronal challenges from a metabolic view. Here we investigated metabolic genes differentially expressed in adult fly glia upon adding neurotransmitters. Searches done using the COG, GO, and KEGG databases identified genes in four types of metabolism including carbohydrates, lipid, amino acid, and alcohol (Table 3). Volcano plot and heatmap were also illustrated in Figure 5, with changes in the expression level of these genes described similarly in Table 1 and Figure 3. Overall, expressions of these metabolic genes were mostly downregulated in ACh- or GABA-treated glia, whereas ATP and Glu also evoked suppression for most genes but with some exception. These results suggest that metabolic suppression might be a general event when glia are in the presence of neurotransmitters. Interestingly, DA triggered increased expression for most metabolic genes, implicating a different response pathway for glia upon DA addition.

Specifically, expressions of the metabolic genes including Hexosaminidase 2 (Hexo2), Trehalose transporter 1-1 (Tret1-1), and Imaginal disc growth factor 2-4, 6 (Idgf2-4, 6) in the Carbohydrate section were downregulated in the ACh- and GABA-treated glia. The expressions of Acyl-CoA synthetase X3 (Acsx3) and a putative lipid phosphatase CG11426 in the Lipid section were downregulated in ACh- and GABA-treated glia, while the expression of CG5254, a putative α-ketoglutarate transporter, was particularly affected by GABA (Figures 5A, D). These results indicate that the DEGs identified for ACh- and GABA-treated glia are similar.

On the other hand, a mixture of upregulated and downregulated genes were identified in the ATP- and Glu-treated glia. For the Carbohydrate section, the expressions of Glycerol-3-phosphate dehydrogenase 3 (Gpdh3), knockdown (kdn), Major Facilitator Superfamily Transporter 12 (MFS12), and SLC22A family member (SLC22A) were upregulated in both types of glia, although in different degrees. The expressions of glucose transporter Glut1 and Insulin-like receptor (InR) were specifically upregulated in Glu-treated glia. The expressions of putative lipid-related genes like CG4586, CG4594, and CG11426 were downregulated in both types of glia, with the expression of CG11426 most affected in the Glu-treated glia. In the Amino acid section, the expressions of Gad1, Eaat2, Glycine transporter (GlyT), and Glycine receptor (Grd) were upregulated in both ATP- and Glu-treated glia (Figures 5B, E).

In the DA-treated glia, most metabolic genes were upregulated in all four types of metabolism groups. The most affected genes included Major Facilitator Superfamily Transporter 9 (MFS9) and SLC22A in the Carbohydrate section; Fatty acid transport protein 3 (Fatp3), pudgy (pdgy), Acyl-CoA synthetase family member 2 (Acsf2), and Acyl-CoA dehydrogenase very long chain (Acadvl) in the Lipid section; Eaat1 and Gs2 in the Amino acid section. Interestingly, sugar transporter 1 (sut1) and the putative L-arginine solute carrier family 25 member 45 (SLC25A45, CG5646) were specifically downregulated in the presence of DA (Figure 5C). Taken together, these results suggest that metabolic suppression in glia is a general event upon neurotransmitter challenges except for DA.

Flies were raised on a fly food (10 g agar, 52.5 g glucose, 60 g dry yeast, 59.1 g corn flour, 26.2 g brown sugar, 10.1 ml antiputrefactive #1 (418 ml propionic acid, 41.5 ml 85% phosphoric acid, and 540.5 ml ddH₂O), and 12.5 ml antiputrefactive #2 (200 g methyl p-hydroxybenzoate add to 95% ethanol to 1 L) in a 12/12 h light/dark condition at 25°C with 70% humidity. All strains were obtained from Bloomington Stock Center, the Vienna Drosophila RNAi Center (VDRC), or as gifts from colleagues. The following stocks were used: repo-GAL4, UAS-mCD8GFP, UAS-mCD8GFP; repo-GAL4, UAS-GCamp6.0, UAS-Gcamp6.0; repo-Gal4.

The flies used to sort glia were kept in incubators at 25°C, 70% humidity, 12/12 h light/dark condition, for 2∼3 days after eclosion. The brains expressing UAS-mCD8GFP; repo-GAL4 were dissected and dissociated as described (Croset et al., 2018; Li et al., 2022; Lu et al., 2023). Male and female adult flies of 100 each were dissected in cold Schneider’s Drosophila medium (SDM: Gibco) and immediately transferred to 1 ml cold SDM with 10 mM of different neurotransmitters added. Neurotransmitters used: acetylcholine chloride (Sigma A6625-25G), dopamine hydrochloride (MCE HY-B0451A), adenosine 5′-triphosphate disodium salt hydrate (Sigma A2383), γ-aminobutyric acid (MCE HY-N0067), and glutamate (Sigma G0355000). After incubated at room temperature on a shaker for 30 min, the brains were washed three times with 1 ml SDM for 10 min each time, then once again with 1XPBS and incubated in 300 μl SDM containing 100 U/ml papain (MCE, HY-P1645) and 2.5 mg/ml collagenase TL (Roche, 5401020001) for 30 min. After incubation for 5 min, the solution was mixed by pipets and centrifuged at 1,500 rpm for another 5 min for three times. A total of 400 μl cold SDM was added to inactivate enzyme. The supernatant was then filtered with 40 μm Cell strainer (Thermo/Fisher, T_70122363547), centrifuged at 400 g, 4°C for 7 min, then re-suspended with 500 μl SDM in a flow tube. Single-cell suspensions were screened by the BD FACS Aria Fusion Cell Sorter flow cytometry and sorted using a 100 μm nozzle. Flow-sorted glia were collected with 500 μl TRIzol Reagent (Thermo, 15596026) in the RNase-free centrifuge tubes. Average 30,000 glial cells were collected in each centrifuge tube and stored in liquid nitrogen.

Total RNA was extracted from the purified glia using TRIZOL and Viral DNA/RNA Mini Purification Kit (MCE, HY-K1082). Then RNA quality was analyzed by the 2100 Bioanalyser (Agilent) and quantified using the ND-2000 (NanoDrop Technologies). Only high-quality RNA sample (OD260/280 = 1.8∼2.2, OD260/230 ≥ 2.0, RIN ≥ 6.5, 28S:18S ≥ 1.0, >1 μg) was used to construct the sequencing library.

RNA-sequencing transcriptome libraries were prepared following the TruSeq™ RNA sample preparation Kit from Illumina using 1 μg of total RNA. Messenger RNA was isolated according to the polyA selection method by oligo(dT) beads and then fragmented by fragmentation buffer. Next, double-stranded cDNAs were synthesized using a SuperScript double-stranded cDNA synthesis kit (Invitrogen, CA) with random hexamer primers (Illumina). The synthesized cDNAs were subjected to end-repair, phosphorylation, and “A” base addition according to the Illumina’s library construction protocol. Libraries were size selected for cDNA target fragments of 300 bps on 2% Low Range Ultra Agarose followed by PCR amplified using the Phusion DNA polymerase (NEB) for 15 PCR cycles. After quantified by TBS380, paired-end RNA-seq sequencing library was sequenced with the Illumina HiSeq xten/NovaSeq 6000 sequencer (2X150 bp read length). The raw paired end reads were trimmed and quality controlled by SeqPrep and Sickle (Pertea et al., 2015). The original data after quality control, namely clean data (reads), were compared with the reference genome to obtain mapped data (reads) for subsequent transcriptome assembly and expression calculation, etc. At the same time, the quality of the comparison results of this transcriptome sequencing was evaluated, software used: HiSat2, TopHat2 (Kim et al., 2015).

In RNA-Seq analysis, gene expression levels were calculated by the number of clean reads (reads counts) located to genomic regions. RSEM software was used to quantitatively analyze the expression levels of genes and transcripts, so as to analyze the differential expression (DE) of genes among different samples (Li and Dewey, 2011). DESeq2 software is used to analyze the DE, according to the comparison to the calculated DEGs read count data. The default screening criteria for significantly DEGs were adjusted using p-value < 0.05 and log2FC ≧ 1, when a gene met these two conditions at the same time, it was considered as a DEG (Love et al., 2014).

With the GO database, genes can be grouped by their function according to the biological processes (BP), cellular component (CC), and molecular function (MF). The DEGs were annotated with GO, and the results were mapped with up-down-regulated gene GO annotation histogram, using software Blast2go (Ding et al., 2018; Liu and Thomas, 2019). Using the KEGG database, genes can be classified according to the pathway they participate. KEGG annotation of DEGs can display the differential genes on the KEGG pathway map, and display the KEGG annotation pathway maps of up-regulated and down-regulated differential genes, using software KOBAS (Kanehisa et al., 2016). Using the EggNOG database, genes can be classified according to their function. COG annotation was made for DEGs displayed on the COG classification statistical map, using the software BLAST+, Diamond (Li et al., 2018).

Quantitative PCR experiments were performed as described (Croset et al., 2018). Total RNAs were extracted from the sorted glia and reverse-transcribed using the HiScript III RT SuperMix for qPCR (+gDNA wiper) (vazyme, R323-01). qPCR was performed in AB Quanut Studio 7 Real-Time PCR System using ChamQ Universal SYBR qPCR Master Mix (vazyme, Q711-02). Student’s t-test was used to statistic differences between control group and treated group, ^nsp > 0.05; *p < 0.05, ^**p < 0.01, ^***p < 0.001, ^****p < 0.0001.

Calcium imaging experiments were performed using adult flies of 3–5 days old (Owald et al., 2015; Jacob and Waddell, 2020). Flies were put on ice for immobilization, then glued into the dish. We used an external saline containing 108 mM NaCl, 5 mM KCl, 2 mM CaCl₂⋅2H₂O, 8 mM MgCl₂⋅6H₂O, 4 mM NaHCO₃, 1 mM NaH₂PO₄, 1 mM NaH₂PO₄⋅2H₂O, 10 mM glucose, 5 mM treharose⋅2H₂O, 10 mM sucrose, 5 mM HEPES, and 1 μM TTX. The fly heads were opened under room temperature in external saline solution. Fluorescence microscope with high frame rate camera was used to record the calcium signals. Images were taken every 40 ms lasting 2 min, and a 10 mM concentration of neurotransmitter was added at 30 s for 5 s. Manual selection and Fiji was used for further analysis (Schindelin et al., 2012). ΔF/F was calculated by fluorescence after addition of neurotransmitter subtracting the baseline fluorescence level (Tian et al., 2009; Vicario and Cali, 2019). To eliminate for baseline differences in F/F0 between samples, we normalized the “preconditioning” treatment for each cell to equal 0 (Park et al., 2022).