Abstract
Homeostatic matching of pre- and postsynaptic function has been observed in many species and neural structures, but whether transcriptional changes contribute to this form of trans-synaptic coordination remains unknown. To identify genes whose expression is altered in presynaptic neurons as a result of perturbing postsynaptic excitability, we applied a transcriptomics-friendly, temperature-inducible Kir2.1-based activity clamp at the first synaptic relay of the Drosophila olfactory system, a central synapse known to exhibit trans-synaptic homeostatic matching. Twelve hours after adult-onset suppression of activity in postsynaptic antennal lobe projection neurons of males and females, we detected changes in the expression of many genes in the third antennal segment, which houses the somata of presynaptic olfactory receptor neurons. These changes affected genes with roles in synaptic vesicle release and synaptic remodeling, including several implicated in homeostatic plasticity at the neuromuscular junction. At 48 h and beyond, the transcriptional landscape tilted toward protein synthesis, folding, and degradation; energy metabolism; and cellular stress defenses, indicating that the system had been pushed to its homeostatic limits. Our analysis suggests that similar homeostatic machinery operates at peripheral and central synapses and identifies many of its components. The presynaptic transcriptional response to genetically targeted postsynaptic perturbations could be exploited for the construction of novel connectivity tracing tools.
SIGNIFICANCE STATEMENT Homeostatic feedback mechanisms adjust intrinsic and synaptic properties of neurons to keep their average activity levels constant. We show that, at a central synapse in the fruit fly brain, these mechanisms include changes in presynaptic gene expression that are instructed by an abrupt loss of postsynaptic excitability. The trans-synaptically regulated genes have roles in synaptic vesicle release and synapse remodeling; protein synthesis, folding, and degradation; and energy metabolism. Our study establishes a role for transcriptional changes in homeostatic synaptic plasticity, points to mechanistic commonalities between peripheral and central synapses, and potentially opens new opportunities for the development of connectivity-based gene expression systems.
- gene expression
- homeostatic plasticity
- synaptic reorganization
- synaptic transmission
- trans-synaptic signaling
- transcriptomics
Introduction
Homeostatic feedback that stabilizes network activity after synaptic weight changes is an important adjunct to correlation-based learning rules (Turrigiano, 2011). Early demonstrations of homeostatic plasticity followed pharmacological manipulations of synaptic transmission in neuronal cultures (Turrigiano et al., 1994, 1998). When global activity levels were artificially increased or decreased, homeostatic forces intervened to maintain firing rates within defined ranges. These homeostatic forces are generated by two processes (Turrigiano, 2011): cell-autonomous changes in intrinsic excitability, which alter the gain of the neuronal voltage response to synaptic currents (Turrigiano et al., 1994; Desai et al., 1999); and adjustments of the synaptic strengths themselves (Petersen et al., 1997; Davis et al., 1998; Turrigiano et al., 1998; Burrone et al., 2002). These adjustments, though in principle achievable in cell-autonomous fashion by altering the density of neurotransmitter receptors in the postsynaptic membrane (Wierenga et al., 2005; Goold and Nicoll, 2010), often involve a trans-synaptic partnership in which postsynaptic neurons communicate deviations from their activity setpoint via retrograde signals to their presynaptic partners, which in turn increase or decrease transmitter release (Cull-Candy et al., 1980; Petersen et al., 1997; Sandrock et al., 1997; Davis et al., 1998; Burrone et al., 2002; Haghighi et al., 2003; Thiagarajan et al., 2005).
Much existing knowledge of retrograde homeostatic communication comes from studies of the neuromuscular junction (NMJ). In mammals and Drosophila, mutations or autoantibodies that reduce the responsiveness of muscle to neurotransmitter cause compensatory increases in motor neuron vesicular release (Cull-Candy et al., 1980; Petersen et al., 1997; Sandrock et al., 1997; Davis et al., 1998). At the Drosophila NMJ, acute pharmacological receptor blockade (Frank et al., 2006) or expression of the inwardly rectifying potassium channel Kir2.1 in muscle (Paradis et al., 2001) induces similar presynaptic compensatory effects. While many gene products and signaling pathways have been implicated in synaptic homeostasis (Davis and Müller, 2015), knowledge of the transcriptional changes that may be required to lock the presynaptic cells into their altered functional state remains scant (Marie et al., 2010).
Pre- and postsynaptic function are also matched at the central synapses between olfactory receptor neurons (ORNs) and projection neurons (PNs) in the antennal lobe of Drosophila (Kazama and Wilson, 2008), where the axons of 20-200 ORNs expressing the same odorant receptor connect to dendrites of an average of three affine PNs in a precise anatomic register (Groschner and Miesenböck, 2019). There is clear covariation between the dendritic arbor sizes of PNs belonging to different transmission channels and the amplitudes of unitary EPSCs: the larger unitary EPSCs of PNs with larger dendritic trees (and, therefore, lower impedances) reflect homeostatic increases in the number of presynaptic ORN release sites in response to increased postsynaptic demand for synaptic drive (Kazama and Wilson, 2008; Mosca and Luo, 2014). This central model of synaptic homeostasis has been characterized physiologically and anatomically, but the molecular mechanism of synaptic matching is unexplored. Taking advantage of the ease with which the presynaptic partners at this synapse can be isolated (they reside in an external appendage, the third antennal segment), we conducted a transcriptome-wide screen for genes regulated by retrograde homeostatic signals. Homeostatic plasticity was induced by adult-onset expression of Kir2.1 in PNs; the expression of a nonconducting mutant of Kir2.1 (Kir2.1-nc) served as control.
Materials and Methods
Drosophila strains and culture
Flies were maintained at 21°C and 65% humidity on a constant 12:12 h light:dark cycle in rich cornmeal and molasses-based food with brewer's yeast. Driver lines GH146-GAL4 (Stocker et al., 1997) and pdf-GAL4 (Renn et al., 1999) were used to target the expression of codon-optimized UAS-Kir2.1 transgenes (see below) to PNs and PDF-expressing clock neurons, respectively. Three copies of two tubulin-GAL80ts insertions on different chromosomes (McGuire et al., 2003) were combined to achieve tight repression of the GAL4-responsive transgenes until induction. The induction incubator was kept at 31°C in 70% humidity on the same 12:12 light:dark schedule.
The cDNA sequence encoding human Kir2.1 was codon-optimized for Drosophila (GenBank accession number MW088713), synthesized at MWG Eurofins, and fused to a codon-optimized N-terminal EGFP tag. The nonconducting variant (Kir2.1-nc) was created by mutating codon 146 of the ion channel sequence from glycine to serine (GGA to AGC) (Haruna et al., 2007). The channel constructs replaced the mCD8::GFP coding sequence in derivatives of plasmid pJFRC2-10XUAS-IVS-mCD8::GFP (Pfeiffer et al., 2010), which were inserted into the attp2 landing site on the third autosome.
Confocal microscopy
Female flies 5 d of age were anesthetized on ice and dissected in PBS (1.86 mm NaH2PO4, 8.41 mm Na2HPO4, 175 mm NaCl). Immediately after dissection, brains were fixed in ice-cold PBS containing 4% (w/v) PFA for 1-2 h at room temperature, rinsed 3 times in ice-cold PBS containing 0.1% (w/v) Triton X-100 (PBT), washed 3 times for 20 min in ice-cold PBT, and mounted and cleared in Vectashield (Vector Labs). Confocal image stacks with an axial spacing of 1-1.5 μm were collected on a Leica Microsystems TCS SP5 microscope with an HCX IRAPO L 25×0.95 W objective.
Electrophysiology
Targeted whole-cell patch-clamp recordings from the fluorescent somata of PNs expressing EGFP::Kir2.1 or EGFP::Kir2.1-nc were obtained through a small cranial window in 5-day-old females. The brain was continuously superfused with extracellular solution containing 103 mm NaCl, 3 mm KCl, 5 mm TES, 8 mm trehalose, 10 mm glucose, 7 mm sucrose, 26 mm NaHCO3, 1 mm NaH2PO4, 1.5 mm CaCl2, 4 mm MgCl2, pH 7.3, and equilibrated with 95% O2–5% CO2. Borosilicate glass electrodes (7-13 MΩ) were filled with intracellular solution containing 140 mm potassium aspartate, 10 mm HEPES, 1 mm KCl, 4 mm Mg-ATP, 0.5 mm Na3GTP, 1 mm EGTA, pH 7.3. Signals were acquired with a MultiClamp 700B Microelectrode Amplifier, filtered at 6-10 kHz, and digitized at 10-20 kHz with an ITC-18 data acquisition board controlled by the Nclamp and NeuroMatic packages. Data were analyzed with NeuroMatic (http://neuromatic.thinkrandom.com) and custom procedures in Igor Pro (WaveMetrics) (Donlea et al., 2014). The membrane time constant was determined by fitting a single exponential to the voltage deflection caused by a 200-ms-long hyperpolarizing current pulse. Input resistances were estimated from linear fits of the subthreshold voltage deflections elicited by 5 pA current pulses of increasing amplitude and a duration of 1 s. Firing rates were quantified by holding cells at resting potentials of –60 ± 2 mV and injecting sequences of depolarizing current pulses (5 pA increments, 1 s duration). Spikes were detected by finding minima in the second derivative of the membrane potential record. The spike rate was calculated by dividing the number of action potentials discharged by the time elapsed between the first and last spike. The current amplitude at which each cell reached a given frequency threshold (1-50 Hz) was used to construct cumulative distribution functions. The distributions were fit with logistic Naka-Rushton functions of the following form (Donlea et al., 2014):
where F is the percentage of cells reaching threshold at a given current level I, Fmax is the percentage of cells reaching threshold at maximal current, I50 indicates the half-maximal or semisaturation current, and the exponent n determines the steepness of the curve. With only two free parameters (I50 and n, given that Fmax is measured experimentally), this simple model provided a satisfying fit to all distributions.
Circadian behavior
Three-day-old female flies were individually inserted into 65 mm glass tubes and loaded into the Trikinetics Drosophila Activity Monitoring system, which was operated at 31°C in 24 h dark conditions for 5–7 d. Group sizes for activity measurements (16 experimental and 16 control flies) reflect the capacity of the monitors.
Third antennal segment dissection
Groups of 20-30 flies were aged in precisely controlled temperature conditions for 5 d (see Fig. 4A) and decapitated with a surgical scalpel on a CO2 pad; the heads were transferred to Petri dishes kept on dry ice. Once a Petri dish contained ∼50 heads, it was sealed with Parafilm and stored at –80°C until RNA extraction. The sealed Petri dishes were dipped in liquid nitrogen for 60 s, vortexed at full strength for 60 s, and then unsealed and placed on a dry-ice-chilled glass stand under a dissection microscope. Individual third antennal segments were picked with fine forceps and placed directly into 100 µl TRIzol (Thermo Fisher Scientific).
RNA extraction
Third antennal segments in 100 µl TRIzol were disrupted with several strokes in a Dounce homogenizer. The homogenates were diluted with 900 µl TRIzol and incubated at room temperature for 5 min. Samples destined for 3′ digital gene expression profiling (3′ DGE) underwent phase separation after the addition of 225 µl chloroform; RNA in the aqueous phase was precipitated with isopropanol and resuspended in 5 µl RNase-free water. Total RNA for RNA-seq and qRT-PCR was isolated with the help of RNeasy minelute columns (QIAGEN), following the addition of 400 µl of 70% RNase-free ethanol to the TRIzol homogenates and on-column DNaseI digests. Samples were snap frozen in liquid nitrogen and stored at –80°C.
cDNA library generation
Libraries for 3′ DGE were generated at MWG Eurofins Genomics from ultrasonically fragmented poly(A)-tailed RNA, which was isolated using oligo(dT) chromatography. Following ligation of an RNA adapter to the 5'-end, the mRNA fragments were reverse-transcribed from an oligo(dT) primer, and the resulting cDNA was PCR-amplified with a high-fidelity polymerase. Each cDNA library was purified, size-selected, quality-checked by capillary electrophoresis, and sequenced on the HiSeq2000 platform (Illumina) in 1 × 100 bp run mode.
For RNA-seq and qRT-PCR, oligo(dT)-enriched RNA underwent 14 cycles of amplification using the SMARTer Ultra Low RNA Kit for Illumina Sequencing (Clontech). After cDNA fragmentation, libraries were prepared in an additional 15 amplification cycles using the NEBNext Ultra DNA Library Prep Kit for Illumina (New England Biolabs) and sequenced on the HiSeq2000 platform (Illumina) in paired-end mode.
Transcriptome analysis
Raw reads were 100 bp in length (paired-end reads for RNA-seq and single-ended reads for 3′ DGE). Fastq files, containing reads and quality scores, were first run through the FastQC package (Andrews, 2010). Highly abundant sequences that did not map to the Drosophila genome (and originated from primers or amplification artifacts) were eliminated using Trimmomatic software (Bolger et al., 2014). Reads were scanned with a 4 bp sliding window and cut when the average quality dropped to <15; trimmed reads <25 bp were discarded. The reads were mapped to Ensembl DM genome release 5.74 using TopHat2 (D. Kim et al., 2013), assigned to transcripts annotated in the transcript file of the Berkeley Drosophila Genome Project release 5.74 with Cufflinks, and merged into an experiment-wide gtf file with Cuffmerge (Trapnell et al., 2012). The gtf file was used to produce raw read counts (using HTSeq) suitable for differential expression analysis in DESeq2 (Love et al., 2014). The topGO and ViSEAGO packages were used to analyze the enrichment of gene ontology (GO) terms in the set of differentially expressed genes called by DESeq2 (unadjusted p < 0.05) vis-à-vis a reference set of all genes with a normalized expression level above 1 (the “gene universe”) (Alexa et al., 2006; Brionne et al., 2019). To keep the number of Fisher's exact tests to a minimum, only GO terms with >40 attached genes were considered. Enriched GO terms with unadjusted p < 0.01 were clustered hierarchically according to Wang's distance, a measure of semantic similarity (J. Z. Wang et al., 2007; Brionne et al., 2019).
qRT-PCR
Transcript levels were determined by qRT-PCR on a LightCycler 480 system (Roche Diagnostics) using SYBR Green I Master Mix (Roche Diagnostics) in 10 µl reactions containing 100 nm of each gene-specific primer and 50 ng of preamplified cDNA. Two sets of primers were designed for each gene of interest. All samples were run in technical triplicates; non–reverse-transcribed mRNA and water served as negative controls. Melting curves were analyzed after amplification, and amplicons were visualized by agarose gel electrophoresis to confirm primer specificity. Relative transcript levels were estimated with the help of the 2–ΔΔCt method (Livak and Schmittgen, 2001), using the housekeeping gene CycK for normalization.
Experimental design and statistical analysis
Sample sizes for differential gene expression analyses were determined in a pilot experiment comparing the expression of circadian-regulated genes during the day and night (Fig. 1D). No statistical tests were used to predetermine sample sizes for electrophysiology, behavior, or qRT-PCR.
Statistical analyses of differential gene expression data and GO term enrichment relied on built-in functions in DESeq2, topGO, and ViSEAGO, as described above (Alexa et al., 2006; Love et al., 2014; Brionne et al., 2019). Hypotheses about differentially expressed genes were evaluated by Wald test in DESeq2 (Love et al., 2014); enrichment of GO terms was quantified using Fisher's exact test in ViSEAGO (Brionne et al., 2019). Group means of electrophysiological parameters, which Shapiro–Wilk tests confirmed were normally distributed, were compared by two-tailed t test.
Transcriptome data are available from the European Bioinformatics Institute's ArrayExpress archive under accession numbers E-MTAB-10 062 (antennal gene expression during the day and night) and E-MTAB-10 065 (antennal gene expression following the induction of Kir2.1 or Kir2.1-nc in PNs).
Results
Antennal transcriptomics
To characterize gene expression in the third antennal segment, 5-day-old male Canton-S flies were decapitated either between zeitgeber time 5 (ZT5) and ZT8 (the day group) or between ZT17 and ZT20 (the night group). After snap-freezing, third antennal segments were manually isolated, and total RNA was extracted in a single batch to minimize variability (Fig. 1A; see Materials and Methods). For both day and night conditions, three biological replicates were prepared, and the resulting six cDNA libraries were sequenced on one lane of an Illumina HiSeq2000 machine using 3′ DGE technology. After stringent quality assessment and read trimming (Fig. 1A), the high-quality reads were mapped to the Drosophila genome (for mapping statistics, see Table 1). Biological replicates showed high correlations with one another (Fig. 1B; Table 1), and day and night samples could easily be distinguished on the basis of their top two principal components (Fig. 1B, inset). Underlying this clean separability were 128 differentially expressed genes, identified by DESeq2 (Love et al., 2014) with a false discovery rate (FDR)-adjusted significance level of < 0.20, and large expression level differences between the day and night (Fig. 1C; Table 2). Core clock components, such as cryptochrome, Clock, period, timeless, and vrille, were found near the top of the amplitude distribution of oscillating transcripts, in two groups at opposite poles of the 24 h cycle (Fig. 1C), consistent with their antagonistic roles in the transcriptional feedback oscillator (Claridge-Chang et al., 2001; McDonald and Rosbash, 2001).
Antennal transcriptomics: workflow, diagnostics, and functional validation. A, Experimental workflow. B, Scatterplot of gene expression levels in biological replicates N1 versus N2. Inset, D and N samples in a principal component analysis (PCA) plot. C, MA plot of log2 fold change in expression (M) versus mean expression level (A) of all transcripts. Triangles represent data points outside the plotted range. D, Expression levels of all transcripts with FDR-adjusted p < 0.20 during the day and night. Each column represents a sequencing library generated from third antennal segments. Red type indicates core clock components.
Mapping metrics of third antennal segment transcriptomes collected during the day (D) and night (N)a
Circadian-regulated genes (FDR-adjusted p < 0.01)a
Transcripts encoding olfactory, gustatory, and ionotropic receptors (ORs, GRs, and IRs) provided an index of the purity of our library preparations. Certain ORs, IRs, and GRs are expressed in antennal ORNs but not elsewhere, while others are absent from antennal ORNs but present in different types of sensory neuron (Clyne et al., 1999; Gao and Chess, 1999; Vosshall et al., 1999; Clyne et al., 2000; Dunipace et al., 2001; Scott et al., 2001; Benton et al., 2009). We detected the former, but not the latter, members of all three receptor families, including the obligatory OR coreceptor Orco/Or83b (Larsson et al., 2004), which was expressed at a level above that of all other OR genes (mean rLog ± SEM = 13.6178 ± 0.1650), in abundance (Fig. 2A-C). These data indicate little, if any, nonantennal contamination of our libraries, and they suggest that receptor gene expression was stable across all experimental conditions.
Expression levels of ORs (A), IRs (B), and GRs (C) in sequencing libraries generated from third antennal segments. Each column represents a library generated during the day (D) or night (N), or after 12, 48, or 96 h induction of Kir2.1 or Kir2.1-nc. Within each induction period, gray dividers separate libraries obtained from flies expressing Kir2.1 (to the left of the divider) from those of flies expressing Kir2.1-nc (to the right of the divider). The gene encoding the obligatory OR coreceptor Orco/Or83b was expressed at a level above those of other OR genes (mean rLog ± SEM = 13.6178 ± 0.16503) and omitted from A.
Transcriptomics-friendly manipulation of postsynaptic excitability
Kir2.1, an inwardly rectifying potassium channel, decreases the input resistance of neurons and clamps their membrane potential at or below resting value; it is widely used as a neuronal “silencer” (Johns et al., 1999; Paradis et al., 2001; Burrone et al., 2002). Some single amino acid substitutions in the P-loop signature sequence of the channel (Heginbotham et al., 1994), such as G146S (here called Kir2.1-nc), block ion flow without affecting the protein's localization (Haruna et al., 2007). We generated Drosophila codon-optimized UAS-EGFP::Kir2.1 and UAS-EGFP::Kir2.1-nc lines and crossed them to the GH146-GAL4 driver, which directs transgene expression to PNs (Stocker et al., 1997). Whereas the nonconducting Kir2.1-nc variant proved innocuous, the expression of functional Kir2.1 under GH146-GAL4 control caused early larval lethality, but this premature death could be circumvented with three tubulin promoter-driven copies of the temperature-sensitive repressor of GAL4, GAL80ts (McGuire et al., 2003), which kept the expression of the channel at bay until the block was thermally relieved during adulthood.
Following their induction for 24 h at 31°C, both EGFP-tagged channels (Kir2.1 and Kir2.1-nc) were detected in PNs of 5-day-old adults at comparable levels and in the same anatomic distribution (Fig. 3A). Whole-cell current-clamp recordings showed that EGFP::Kir2.1 lowers the input resistance and membrane time constant relative to EGFP::Kir2.1-nc (Fig. 3B-D) and powerfully opposes depolarization: Kir2.1-expressing neurons required approximately twofold larger depolarizing currents to drive spiking across a firing rate range of 1-50 Hz (Fig. 3E). Although Kir2.1 does not strictly silence the population of neurons in which it is expressed (the added potassium conductance can always be compensated by a large enough current injection; Fig. 3B,E), the currents necessary to do so seem difficult to attain in vivo.
A transcriptomics-friendly neuronal activity clamp. A, Maximum intensity projections of confocal image stacks through the antennal lobes of 5-day-old female flies carrying EGFP::Kir2.1 or EGFP::Kir2.1-nc transgenes under GH146-GAL4 and tub-GAL80ts control. The expression of Kir2.1 constructs is undetectable at 21°C (top) but induced at 31°C (center and bottom). Scale bar, 20 µm. B, Example voltage responses to 5 pA current steps of antennal lobe PNs expressing EGFP::Kir2.1-nc (black) or EGFP::Kir2.1 (red). C, D, Kir2.1 (red) lowers the input resistance Rm (t(14) = 2.1652, p = 0.0481; C) and shortens the membrane time constant τm (t(14) = 4.4959, p = 0.0005; D) relative to Kir2.1-nc (black). Circles represent individual PNs. Error bars indicate mean ± SEM. E, Cumulative distribution functions of the percentages of PNs reaching a spike frequency of 30 Hz at different levels of injected current (left); semisaturation currents (middle) and percentages of cells reaching spike rates of 1-50 Hz, for PNs expressing Kir2.1-nc (black) or Kir2.1 (red). F, Circadian locomotor rhythms in constant darkness. Locomotion was quantified as the total number of midline crossings per minute in groups of 16 flies expressing Kir2.1-nc (black) or Kir2.1 (red) under pdf-GAL4 control. The traces were smoothed with a Gaussian kernel (1.25 h FWHM) and show data collected on days 2-5 after the flies were transferred to activity monitors.
A simple behavioral test supported this conclusion. Adult-onset expression of Kir2.1 in the PDF-expressing ventral subset of lateral pacemaker neurons (using the pdf-GAL4 driver) (Renn et al., 1999) disrupted the circadian locomotor rhythm in constant darkness, as expected (Nitabach et al., 2002), whereas flies expressing Kir2.1-nc remained rhythmic (Fig. 3F).
Trans-synaptic regulation of gene expression: transmitter release and synapse remodeling, and a late shift to proteostasis and neuroprotection
To delineate changes in presynaptic gene expression after muting postsynaptic neural activity, we compared the third antennal segment transcriptomes of flies expressing either Kir2.1 or Kir2.1-nc in PNs (Fig. 4A). We studied three induction times (12, 48, and 96 h) in individuals that were age-matched at the point of analysis: all tissues were harvested between ZT6 and ZT7 on the fifth posteclosion day (Fig. 4A). Two sequencing technologies (3′ DGE for the 12 and 48 h groups and standard RNA-seq for the 96 h group) gave similar mapping metrics (Tables 3 and 4).
Trans-synaptic regulation of gene expression: transmitter release and synapse remodeling, and a late shift to proteostasis and neuroprotection. A, Experimental design. B, Enrichment of GO biological process terms in third antennal segment transcriptomes after 12, 48, and 96 h induction of Kir2.1. The dendrogram represents semantic groupings among GO terms. C, Expression levels of transcripts attached to two semantic groupings, “transmitter release” and “synapse remodeling” (B), after 12 h of induction of Kir2.1 or Kir2.1-nc. Each column represents a sequencing library. Red type indicates gene products previously implicated in homeostatic synaptic plasticity (see Discussion).
Mapping metrics for third antennal segment transcriptomes collected after 12 or 48 h of induction of Kir2.1 (K) or a nonconducting control (C)a
Mapping metrics for third antennal segment transcriptomes collected after 96 h of induction of Kir2.1 (K) or a nonconducting control (C)a
For 12 h induction, experimental (Kir2.1) and control (Kir2.1-nc) flies were placed at 31°C from ZT18 until ZT6 on their fifth posteclosion day and decapitated between ZT6 and ZT7 on the same day (Fig. 4A). Three biological replicates were sequenced for each genotype: one from males and two from females. The inclusion of sex as a variable allowed us to verify that modest expression level changes of low-abundance transcripts could be detected (Shiao et al., 2013). A DESeq2 differential expression analysis based on sex alone returned many known sexually dimorphic ORs at an FDR-adjusted p < 0.20, including Or47b (upregulated in males) (Kopp et al., 2008; Shiao et al., 2013), Or85a (upregulated in females) (Kopp et al., 2008), and Or22b (upregulated in females) (Kopp et al., 2008). After sex differences were accounted for and removed by the regression model entered into DESeq2, a comparison of flies expressing Kir2.1 versus Kir2.1-nc in PNs highlighted 25 differentially expressed third antennal segment genes with FDR-adjusted p < 0.20 (Table 5). The average changes in absolute expression levels of the top 20 differentially expressed genes were ∼15-fold smaller than those of the top 20 clock-controlled genes (log2 fold changes: 0.35 ± 0.08 for homeostatic genes vs 1.9 ± 0.9 for circadian genes; Tables 2 and 5), resulting in many fewer significant hits for the same FDR threshold. Among genes with the smallest FDR-adjusted p values, many are involved in cell fate commitment and morphogenesis (Table 5); eight (bazooka, sugar-free frosting, plum, prospero, Ankyrin 2, spätzle, Syncrip, and ATP6AP2) have been linked to synaptic organization or synapse formation, however indirectly (Doe et al., 1991; Ruiz-Canada et al., 2004; Koch et al., 2008; Pielage et al., 2008; Baas et al., 2011; Sutcliffe et al., 2013; Yu et al., 2013; Halstead et al., 2014; Dubos et al., 2015).
Differentially expressed genes after 12 h of Kir2.1 induction (FDR-adjusted p < 0.20)a
For 48 h induction, experimental and control groups were shifted to 31°C at ZT6 of their third posteclosion day and decapitated between ZT6 and ZT7 on day 5 (Fig. 4A). Three biological replicates (all from males) were sequenced for each genotype using 3′ DGE technology. One of the Kir2.1-nc replicates (C48-2) did not cluster well with the others (Table 3) and was excluded from the differential expression analysis, which produced 26 hits with FDR-adjusted p < 0.20 (Table 6). Conspicuous among these hits were several ribosomal components and three chaperones of the Hsp20 family (Hsp27, Hsp67Bc, and Hsp23) (Haslbeck et al., 2019). At first glance, the upregulation of heat shock proteins might suggest a direct effect of our method of transgene induction (31°C heat), but on reflection heat cannot explain the observed differences because experimental and control flies were exposed to the same temperature regimen. A more plausible explanation is, therefore, that prolonged postsynaptic silencing places an intense homeostatic burden on presynaptic partners which elicits a generalized increase in protein synthesis.
Differentially expressed genes after 48 h of Kir2.1 induction (FDR-adjusted p < 0.20)a
For 96 h induction, experimental and control groups were kept at 31°C from ZT6 of their first posteclosion day and again decapitated between ZT6 and ZT7 on day 5 (Fig. 4A). A total of 20 libraries were sequenced in two batches using RNA-seq technology. A different sequencing method was chosen to ensure that our results were valid across sequencing platforms, and more replicates were processed to increase sensitivity. The first batch consisted of 12 samples with six replicates from each of the two genotypes (all male third antennal segments). Two replicates of each genotype in the first batch (K96-2, K96-3, C96-2, and C96-3) were sequenced to twice the depth of the others to detect very lowly expressed genes more reliably. The second batch (eight samples in total) consisted of another four samples of each genotype: two from females and two from males. Two samples (K96-1 and K96-9) had low within-batch correlations and were omitted from the analysis (Table 4). The increase in statistical power enabled the detection of 32 differentially expressed genes with FDR-adjusted p < 0.05 after controlling for sex and batch in DESeq2 (Table 7). Three biological processes stand out among these differentially expressed genes. First, six genes related to the Imd and Toll pathways of the innate immune response (Valanne et al., 2011) were strongly downregulated: the pattern recognition receptor PGRP-SD; the antibacterial peptide Drosocin (Dro); the negative regulator of Imd, pirk; and the antimicrobial peptides Bomanin Short 1, 3, and 5 (also known as IM1, IM2, and IM3). Second, chaperones of the Hsp20 family, already encountered after 48 h induction, were again upregulated (Hsp26, Hsp23, and Hsp67Bc) (Haslbeck et al., 2019). And third, four genes involved in programmed cell death were differentially expressed, with two pro-apoptotic factors downregulated—matrix metalloproteinase 1 (Mmp 1) and apoptosis-inducing factor (AIF)—and two gene products inhibiting apoptosis upregulated (Hsp26, Buffy) (Quinn et al., 2003; H. D. Wang et al., 2004; Joza et al., 2008). Overall, the 96 h picture suggests a transcriptional landscape skewed toward cell protection and maintenance.
Differentially expressed genes after 96 h of Kir2.1 induction (FDR-adjusted p < 0.05)a
Enriched GO biological process terms after 12 h of Kir2.1 inductiona
To obtain an aerial view of transcriptionally regulated biological processes during all induction periods, we probed for coordinated changes in functionally related sets of genes via GO enrichment analyses. These analyses were performed on all differentially expressed genes with unadjusted p < 0.05 and included only GO terms with >40 attached genes; the enriched GO terms (Fisher's exact test, p < 0.01) were then hierarchically clustered according to their semantic similarity (J. Z. Wang et al., 2007; Brionne et al., 2019) (Fig. 4B). After 12 h of induction, presynaptic transcriptional changes centered on genes encoding synaptic release and remodeling machinery; at 48 h and beyond, protein synthesis and degradation, and energy metabolism, predominated (Fig. 4B,C; Tables 8–10). Closer scrutiny of the 81 genes responsible for the early enrichment of synaptic GO annotations (Fig. 4C; Tables 11 and 12) uncovered many with established roles in homeostatic plasticity at the NMJ (or with known interactions with such genes), as we discuss below. Although typical transcripts showed only modest expression level changes of 15%-30%, their regulation was clearly visible across multiple libraries (Fig. 4C). This consistency across biological replicates, and the statistically verified overabundance of synaptic genes in the differentially expressed set with low unadjusted p values (Tables 11 and 12), suggest a genuine signal.
Cross-validation of regulated genes with 3′ DGE, RNA-seq, and qRT-PCR
As a further validation of our gene expression measurements, we compared transcriptome-wide 3′ DGE with transcriptome-wide RNA-seq data. There was an approximately linear relationship between the average expression levels of all genes in all samples (Fig. 5A), with a small departure in lowly expressed genes caused by the extra amplification step in the RNA-seq protocol (see Materials and Methods); as a result, RNA-seq reported systematically higher expression levels for scarce transcripts than did 3′ DGE. For genes transcribed at moderate to high expression levels, the two sequencing platforms were in close agreement.
Cross-validation of differential gene expression. A, Scatterplot of average gene expression levels determined by RNA-seq versus 3′ DGE. B, Scatterplot of log2 fold changes (mean ± SEM) in the expression levels of 11 transcripts after 96 h induction of Kir2.1, determined by RNA-seq or qRT-PCR.
We next selected 11 transcripts for qRT-PCR verification. These transcripts were chosen from the set of differentially expressed genes (unadjusted p < 0.05; both upregulated and downregulated) in the 96 h induction group and used to validate all deep sequencing data. The fold changes of the 11 chosen transcripts, as estimated by qRT-PCR with normalization to the housekeeping gene CycK, correlated tightly with 3′ DGE and RNA-seq measurements (Fig. 5B). This agreement between three independent measures of gene expression, at a transcriptome-wide scale and across several individual genes, lends confidence to our analysis.
Discussion
Trans-synaptic regulation of gene expression
Our study introduces an experimental system for detecting changes in gene expression in response to changes in the electrical excitability of a partner cell. The product of the Kir2.1 transgene powerfully suppresses the activity of neurons in which it is expressed, while a control transgene, which codes for a potassium channel with a single amino acid substitution in its selectivity filter, has no effect (Fig. 3B-F). Isogenic strains expressing one or the other of these transgenes from the same chromosomal locus offer an ideal platform for differential gene expression analyses because differences between them can be pinned to a single codon change in the genome. The finding that prolonged postsynaptic silencing induces the expression of Hsp20 proteins in a manner unrelated to heat shock (Tables 6 and 7) underscores the power of this carefully controlled system.
The same finding, however, also highlights a limitation particular to our current approach. We imposed the Kir2.1 clamp on the first synaptic relay in the Drosophila olfactory system because its presynaptic and postsynaptic elements are easily separable by purely physical means, but this convenience exacted a price: the third antennal segment contains not only ORNs but also glial and support cells, which account for approximately two-thirds of the segment's cell population (Vosshall et al., 1999). We are therefore unable to determine whether the expression of Hsp20 proteins is exclusively or even partially neuronal. Although the same reservation does not apply to the many synaptic genes that are differentially expressed during the early phase of the homeostatic response (Fig. 4B,C), the presence of non-neuronal elements may nevertheless have hindered the detection of low-abundance neuronal transcripts or underestimated their fold change. Both of these drawbacks could be overcome by FACS isolation of a genetically labeled cell population from dissociated neural tissue, as would be required as a matter of course in all instances where the synaptic partners are intermingled. With this extra step (which we were unable to take because antennal tissue could not be recovered intact from its chitinous shell), our system will be easily adapted for analyses of transcriptional changes elicited in presynaptic cells by a loss of postsynaptic responsivity, in postsynaptic cells by a loss of presynaptic input, or in glial cells by a heightened demand for synaptic remodeling.
Despite these caveats, many of the early expression level changes we detect affect genes encoding synaptic proteins with known, suspected, or at least plausible roles in homeostatic plasticity (Fig. 4B,C; Tables 5, 11, and 12) (Davis and Müller, 2015): elements of the wingless signaling system (e.g., Wnk, sgg), which acts as an endogenous suppressor of homeostatic compensation at the NMJ (Marie et al., 2010); the v-SNARE synaptobrevin (nSyb) and its chaperone Nsf2 (Söllner et al., 1993; Bacci et al., 2001); rab3 guanine nucleotide exchange factor (rab3-GEF), which controls the assembly and distribution of active zone components (Bae et al., 2016) and regulates the nucleotide state-dependent association of rab3 with synaptic vesicles, which in turn determines the calcium sensitivity of their release (Geppert et al., 1997; Müller et al., 2011); an active zone resident (unc-13-4a) known to associate with the Rab3-interacting molecule RIM and other active zone components (Schoch et al., 2002; Liu et al., 2011; Müller et al., 2012); a kinesin motor heavy chain (Khc-73) implicated in active zone assembly and synaptic homeostasis (Tsurudome et al., 2010); an active zone-integral guanylate kinase (CASK) that serves as a phosphorylation target of CDK5 (Samuels et al., 2007), which homeostatically regulates presynaptic calcium influx and release probability (Seeburg et al., 2008; S. H. Kim and Ryan, 2010); the E3 ubiquitin-protein ligase highwire (hiw) and the Smad protein Medea (Med), which in motor neuron terminals are part of the transduction cascade for a retrograde signal from muscle (Haghighi et al., 2003; McCabe et al., 2004; Goold and Davis, 2007); the cytoskeletal anchor Ankyrin 2 (Koch et al., 2008; Pielage et al., 2008); and subunits or accessory proteins of voltage-gated ion channels (quiver, ether-á-go-go, Hyperkinetic, paralytic) (Tables 5, 11, and 12). Collectively, these changes could signal an increase in the number of release sites or an expansion of the release-ready vesicle pool, inferred to represent the dominant quantal parameter change during homeostatic matching at ORN-to-PN synapses (Kazama and Wilson, 2008) and one of two homeostatic levers at the NMJ (the other being modulation of calcium influx into the terminal) (Müller et al., 2012).
When drawing comparisons with earlier work, however, it is important to bear in mind differences in the speed of induction and expression of the homeostatic response. Abrupt adult-onset PN silencing resembles an acute postsynaptic receptor blockade at the NMJ more closely than it does the slow developmental processes studied in analyses of arbor size matching in the antennal lobe (Kazama and Wilson, 2008; Mosca and Luo, 2014), but homeostatic compensation at the NMJ is evident within minutes, long before changes in gene expression can occur (Frank et al., 2006). That elements of the homeostatic machinery are encoded by trans-synaptically regulated genes must therefore reflect a secondary layer of feedback or the more profound reallocation of ORN synapses between PNs and other postsynaptic partners, such as local neurons of the antennal lobe (Groschner and Miesenböck, 2019).
Because changes in the expression levels of putative homeostatic genes are small compared with those of circadian-regulated genes (Figs. 1D, 4C), we were forced to relax our FDR thresholds in the 12 and 48 h induction experiments, raising the specter of false positives in these datasets. Two observations should allay this concern. First, the 96 h induction experiment, whose greater statistical power allowed a more stringent significance threshold, recovered many of the same biological processes and indeed the same genes (e.g., Hsp23, Hsp67Bc) as the statistically weaker 48 h induction experiment (Tables 9 and 10). Second, our qRT-PCR validation included several genes that failed to cross the most stringent FDR threshold (Fig. 5B). These qRT-PCR spot checks confirmed that expression level changes detected by RNA-seq or 3′ DGE were accurate. Nonetheless, new candidates emerging from our screen will need to survive rigorous functional studies before joining the ranks of established homeostatic plasticity genes.
Enriched GO biological process terms after 48 h of Kir2.1 inductiona
Enriched GO biological process terms after 96 h of Kir2.1 inductiona
Labeling connections with trans-synaptically regulated genes?
In many neurobiological studies, the object of interest is not a population of genetically homogeneous neurons but an operational unit (a circuit) defined by connectivity rather than a common genetic marker. Circuit analyses have benefited greatly from the development of trans-synaptic vectors, which travel along synaptic connections between specific types of neuron and serve as vehicles for the distribution of other encodable tools (Miesenböck and Kevrekidis, 2005; Miesenböck, 2009; Sjulson et al., 2016). Ideally, trans-synaptic expression systems possess a mechanism that allows their initialization at a specific location, a rule that governs their propagation in the network, and gain. Viruses have some of these characteristics (Strack et al., 1989; Ugolini et al., 1989; Card et al., 1990; Wickersham et al., 2007). Their infectious spread can follow routes of synaptic transmission, and replicative gain (where permitted) allows each infected neuron to supply more viral particles to its outputs than it receives from its inputs. Viral infections are, however, difficult to control and initialize with single-cell resolution and can produce considerable toxicity and extrasynaptic spread. Clearly, the ideal trans-synaptic vector would, instead of carrying its own genetic material, act on expression cassettes that lie dormant in the genome of the host organism until switched on by a trans-synaptic signal.
Circuit-tracing systems such as trans-Tango, TRACT, and BAcTrace are built on this principle but require the reconstitution of an exogenous cell-to-cell signaling apparatus (Huang et al., 2017; Talay et al., 2017; Cachero et al., 2020). This adds genetic complexity and the danger of overexpression artifacts if the foreign molecules escape synaptic confinement. Eavesdropping on endogenous trans-synaptic communication during homeostatic plasticity offers a possible cure for these problems. Imagine a sudden, targeted loss of excitability in a small group of neurons or even a single cell, brought about by the inducible expression of Kir2.1. If presynaptic partners compensate for this perturbation, the upregulation of plasticity genes could be coupled to the expression of sensors, actuators, transcription factors, or recombinases.
The chief obstacle to the development of this retrograde tracing technology is the small, at most twofold, changes in homeostatic gene expression we detect (Tables 5, 11, and 12). We suspect that these changes will need to be amplified with adequate signal-to-noise ratio, perhaps by flipping a recombination switch, to be practically useful. The application of this strategy to immediate early genes serves as an encouraging precedent: immediate early genes are used widely to trap neural ensembles in defined functional contexts, although their activity-dependent expression level changes roughly equal those of our trans-synaptically regulated genes (Sjulson et al., 2016; DeNardo and Luo, 2017). And while the degree or direction of homeostatic adjustment may vary among synapses (Turrigiano, 2011), the substantial overlap between elements of the homeostatic machinery at the NMJ (Davis and Müller, 2015) and homeostatically regulated genes in the antennal lobe (Fig. 4C) points to a conserved mechanistic core.
Differentially expressed genes (p < 0.05) attached to GO biological process terms in the semantic grouping “transmitter release” after 12 h of Kir2.1 induction
Differentially expressed genes (p < 0.05) attached to GO biological process terms in the semantic grouping “synapse remodeling” after 12 h of Kir2.1 induction
Footnotes
This work was supported by Wellcome grants 209235/Z/17/Z, 106988/Z/15/Z, 090309/Z/09/Z, and 089270/Z/09/Z; Gatsby Charitable Foundation grant GAT3237; and European Research Council grant 832467. We thank Paul Overton for providing advice on RNA isolation; Amélie Baud for giving many helpful tips for the gene ontology analysis; Ruth Brain for assisting with stock maintenance and dissections; Jessica Beevers for helping with dissections and making delicious fly food (from the flies' perspective); and Ronald Davis, Jeffrey Hall, Mani Ramaswami, Gerald Rubin, and Reinhard Stocker for kindly sharing reagents.
The authors declare no competing financial interests.
- Correspondence should be addressed to at Evan R. Harrell at evan-richard.harrell{at}pasteur.fr or Gero Miesenböck at gero.miesenboeck{at}cncb.ox.ac.uk
This is an open-access article distributed under the terms of the Creative Commons Attribution 4.0 International license, which permits unrestricted use, distribution and reproduction in any medium provided that the original work is properly attributed.
