Transcriptome Analysis Reveals the Essential Role of the Domesticated Gag Gene of Drosophila LTR Retrotransposons, Gagr, in the Modulation of Immune Response

1. Introduction

A major factor in the formation of new genes can be the molecular domestication of retroelement sequences, such as retrotransposons and retroviruses. Mammals' domesticated genes of retroelements have been extensively researched; members of various gene families of retroviral origin, such as PNMA, Mart, and syncytins, have been identified and are critical to placenta formation and embryogenesis [1,2]. A number of other genes may be involved in protection against retroviruses entering the body. For example, mice that have 4 such genes – Fv1, Fv4, Rmcf1, and Rmcf2 [3]. Many vertebrate transcription factors have the SCAN domain, which was produced by domesticated sequences of retroelements. [4]. Some domesticated retroelement sequences have recently been characterized: for example, the sequence upstream of the Pparg gene, required for adipogenesis, is derived from the LINE retroelement and plays the role of an enhancer [6]; the PRLH1 transcript originated from the endogenous retrovirus ERV-9 and is involved in the repair of double-strand breaks [6].

The Gagr gene is an illustration of the molecular domestication of the gag gene of retrotransposons/retroviruses in invertebrates. [7]. It has been shown that the role of the Gagr gene may be associated with the immune response and participation in processes that are associated with stress reactions [8]. However, the function of the Gagr gene still remains unknown.

Orthologs of this gene are found in all sequenced genomes of the Drosophila genus; they possess a highly conservative structure and are the result of long-term domestication [7]. There are several research results for the Gagr gene that indirectly indicate its involvement in a number of important processes related to stress reactions. Gagr expression is activated in response to the induction by bacterial lipopolysaccharides in S2 cells, and this activation depends on the regulators of the MAPK/JNK stress signaling pathways Tak1, hep and bsk [9]. Gagr expression increases significantly after intraabdominal injection of DCV viruses (Drosophila C virus), FHV (Flock House virus) and SINV (Sindbis virus) [10].

Gagr protein interactions with other proteins that were established during the identification of protein–protein complexes in S2R+ D. melanogaster cells [11]. The Gagr protein is transmembrane [7] and physically interacts with five partner proteins (14-3-3epsilon, Pdi, eIF3j, CG3687 and CG6013). At least three Gagr partners, 14-3-3epsilon, Pdi and eIF3j, are involved in stress-related functions. 14-3-3epsilon is a conservative regulator of the activity of MAPK and other stress signaling pathways in animals [11]. Chaperone Pdi plays an important role in endoplasmic reticulum stress (ER-stress) and UPR (unfolded protein response) [13]. eIF3j is a subunit of translation initiation factor. eIF3j has been shown to be necessary for IRES-dependent translation which occurs under conditions of cell stress [14]. The CG3687 gene is a poorly studied in D. melanogaster. The only known fact is knockdown of the CG3687 gene results in a flightless phenotype [15]. The last Gagr partner is CG6013 protein with an unknown function, which has homology with the human CCDC124 protein. In the yeast Saccharomyces pombe, the product of the orthologous Oxs1 gene is a cofactor of the transcription factor Pap1 (AP-1-like transcription factor) in the Pap1/Oxs1 signaling pathway [16]. Consequently, it is important to investigate the function of the Gagr gene in cell stress, considering the existing information regarding its activation and protein–protein interactions.

The Gagr gene promoter contains one binding motif for the kayak transcription factor, which is one of the components of the JNK signaling pathway, and two binding motifs for the Stat92E transcription factor, which is a component of the Jak-STAT signaling pathway [8].

One of the methods for studying the function of the Gagr gene is its inactivation. Knockdown allows to specifically turn off gene expression through RNA interference. In the presented work, flies with Gagr gene knockdown in all tissues were obtained. The study examined the lifespan, imago mobility, resistance to heat stress, and the effect of ammonium persulfate (APS), a stressor that disrupts cellular homeostasis, on the transcriptomes of male and female flies with the Gagr gene knockdown.

2. Results

2.1. Physiological tests of the Gagr gene knockdown mutant

2.1.1. Knockdown of the Gagr gene does not affect embryonic and larval mortality of flies

The Tub-GAL4 driver strain is heterozygous for the dominant Tabby mutation, which is localized on chromosome 3; its phenotypic manifestation is a short body. On the same chromosome, in a trans position relative to the mutant Tabby allele, the GAL4 gene is localized, which is necessary for the induction of knockdown. The genetic construct used to drive Gagr gene knockdown is located on the second chromosome of the responder strain. In the crossing of the driver and responder strains, flies with a long body with a knockdown and flies with a short body, as a by-product of crossing, should be obtained in a 1:1 ratio. To test the hypothesis (1:1), we used the χ² method (Table 1). For the control strain, χ²=1.57+1.57=3.14<3.841. In the Gagr knockdown strain, χ²=1.76+1.76=3.52<3.841. Thus, knockdown of the Gagr gene does not affect embryonic and larval mortality of flies.

2.1.2. Knockdown of the Gagr gene affects the lifespan of flies under standard and stress conditions

We measured the lifespan of the Gagr knockdown mutants in comparison to control flies (Figure 1a). The maximum lifespan of males with the Gagr gene knockdown was 45 days. At the same time, the maximum lifespan of Tub-GAL4>w¹¹¹⁸ males was 55 days. The maximum lifespan of Tub-GAL4>UAS-Gagr females was 55 days. At the same time, the maximum lifespan of Tub-GAL4>w¹¹¹⁸ females was 75 days. As a result, Tub-GAL4>UAS-Gagr flies have a shorter maximum lifespan than the control strain.

Next, we measured the survival rate of flies on a medium with 0.1 M APS. The maximum lifespan of males with the Gagr gene knockdown was 30 hours, females – 50 hours (Figure 1b). For Tub-GAL4>w¹¹¹⁸ males, the maximum lifespan was about 50 hours, for females – 64 hours. Thus, Tub-GAL4>UAS-Gagr flies die faster under APS stress conditions than Tub-GAL4>w¹¹¹⁸ flies (both females and males). When exposed to APS, females are more resistant to stress than males.

2.1.3. Knockdown of the Gagr gene does not lead to changes in adult motility

According to Flybase, knockdown of one of the genes, CG3687, whose product is a partner of the Gagr protein, leads to the flightless phenotype (https://flybase.org/reports/FBgn0034097). We supposed that knockdown of the Gagr gene could affect the function of C3687 and somehow affect fly motility. No significant difference in vertical ascent time was found between flies with the Gagr gene knockdown and control flies. All individuals covered a distance of 17 cm in 10 ± 1 sec. Thus, knockdown of the Gagr gene does not lead to a decrease in motility.

2.1.4. Knockdown of the Gagr gene leads to increased resistance to heat stress

The heat stress survival study was performed in a water bath at 38°C (Figure 1c). The best survival was observed in flies with knockdown of the Gagr gene. After 2 hours of heat stress followed by 24 hours of rest, the number of surviving flies was 61.70(±23.86)% for Tub-GAL4>UAS-Gagr hybrids and 35.81(±7.59)% for Tub-GAL4> hybrids w¹¹¹⁸ . Thus, resistance to heat shock correlates with knockdown of the Gagr gene.

2.1.5. Knockdown of the Gagr gene in females promotes the occurrence of melanized masses in the fat body

Despite the absence of morphological manifestations of the Gagr gene knockdown, we found that females had black spots on the body, visible through the cuticle of the abdomen. When flies were dissected, multiple black granules, clusters of melanin capsules, were found in the fat body (Figure 1d). We did not find similar manifestations in the mutant males.

2.2. Transcriptomic analysis of the Gagr gene knockdown mutant

2.2.1. Analysis of differentially expressed genes in the Gagr knockdown strain

We analyzed the transcriptomes of males and females obtained from the Tub-GAL4>w¹¹¹⁸ and Tub-GAL4>UAS-Gagr hybrid flies (Table S1; Figure 2). The transcription level of 297 genes was increased in females with the Gagr knockdown in comparison to control females, and transcription level of 191 was increased genes in males (|Log₂FoldChange|>0.6, Padj<0.05).

Enrichment of genes by functional categories (GO) showed that Tub-GAL4>UAS-Gagr females had significantly increased transcription of genes that respond to infection with Gram-positive and Gram-negative bacteria. In males there are fewer of these categories, but the transcription of immune response genes is also increased.

Transcription of 34 genes was increased in both males and females. Among them there are the genes for antimicrobial peptides: AttA, AttC, AttD, DptB, LysB, and Dro. This means that knockdown of the Gagr gene leads to dysregulation of the Imd and Toll signaling pathways.

The transcription of 263 genes is reduced in Tub-GAL4>UAS-Gagr females relative to control females, and the expression of 153 genes is reduced in Tub-GAL4>UAS-Gagr males. 28 genes have reduced expression in both males and females. Common genes included Lsp1beta, MtnE, Cyp4d1, Yp1, Mal-B1. Enrichment analyses showed that females had reduced expression of detoxification genes and a number of metabolic processes. In males, downregulated genes were associated with digestion, response to zinc, mRNA splicing, and fertilization.

One of the well-studied targets of the Jak/STAT pathway is the TotA gene. Its transcription is activated by the transcription factor STAT92E. This gene had reduced expression in the Tub-GAL4>UAS-Gagr males, and did not differ in expression in females. The pattern of its expression may indicate a reduced level of Jak/STAT pathway activation in males with the Gagr gene knockdown.

Thus, in both males and females with knockdown of the Gagr gene, transcription of the target genes of the Imd and Toll pathways, which respond to Gram-positive and Gram-negative bacteria, was activated, but not the target genes of the Jnk and Jak/STAT pathways. This means that the Gagr gene, which contains sites in the promoter region for binding the transcription factors STAT92E and kayak as a target of the JNK and Jak/STAT pathways, influences the regulation of the Imd and Toll immune response pathways.

Next, we confirmed the transcription level of individual immune response genes using real-time PCR (Figure 3).

The analysis included differentially expressed genes in Tub-GAL4>UAS-Gagr strain: the Gagr gene; the genes of antimicrobial peptides, TotA, TotC, AttB, and CecA2, expressed at a high level in the fat body of flies; the gene of the JNK pathway repressor, Socs16D; the genes of serine endopeptidase inhibitors involved in regeneration processes, Spn88Eb and Spn28Dc (inhibitor of the melanization process); the serine endopeptidase genes, Ser6, CG1304, CG10051 and CG10232; the antiviral effector gene Nazo expressed downstream of Sting and Rel (Relish) signaling, transcription factor of the Imd and Toll pathways; the G33346 gene, predicted to code a protein with RNA and single-stranded DNA endonuclease activity; the Ppo1 gene of the prophenoloxidase involved in the process of melanization.

The studied genes were expressed differently in males and females. The genes of AMPs TotA, TotC, AttB were expressed at a low level in Tub-GAL4>UAS-Gagr males, and CecA2, on the contrary, at a high level. The serpin genes, Spn88Eb and Spn28Dc, were expressed at a higher level in Tub-GAL4>UAS-Gagr males, as well as the genes CG10051, CG10232, GG1304, Ser6. Surprisingly, the Ppo1 gene was expressed at a low level in females and at a high level in males as well as its inhibitor Spn28Dc. The Nazo gene had the highest expression level in mutant females. This data indicates the activation of the Imd and Toll pathways in flies with knockdown of the Gagr gene.

Thus, we confirmed that knockdown of the Gagr gene affects the transcription of immune response genes in males and females.

2.2.2. Analysis of the transcriptomic response to the action of ammonium persulfate

APS exposure resulted in changes in the transcription level of a significant number of genes in females and males (Figure 4). In Tub-GAL4>w¹¹¹⁸ females, 418 genes increased expression and 795 genes decreased expression (|Log₂FoldChange|>0.6, Padj<0.05).

Functional analysis of DEGs showed that the effect of APS is primarily associated with an increase in the expression of genes associated with stress: activation of the immune response (inflammatory response), response to injury, response to stress, as well as genes associated with the regulation of biological processes. The activated genes included mainly components of the Imd and Toll immune signaling cascades. The role of these genes in the response to abiotic stress factors is poorly studied (in contrast to their role in the immune response), however, they are of great interest for identifying the universal mechanisms of the stress response in Drosophila. Genes that decreased their expression in Tub-GAL4>w¹¹¹⁸ females were associated with cellular processes, detoxification and localization.

In Tub-GAL4>w¹¹¹⁸ males, 171 genes increased their expression. Enrichment of upregulated genes also revealed categories associated with immune response. 44 genes decreased their expression in Tub-GAL4>w¹¹¹⁸ males. Enrichment of these genes yielded only one category: metabolic processes.

In Tub-GAL4>UAS-Gagr females, 411 genes increased their expression, 458 genes decreased their expression. Enrichment revealed approximately the same categories as for Tub-GAL4>w¹¹¹⁸ females. In Tub-GAL4>UAS-Gagr males, 67 genes increased their expression and 42 genes decreased their expression. Enrichment of genes revealed categories associated with metabolic processes. If we compare the sets of genes that changed their expression Tub-GAL4>UAS-Gagr and Tub-GAL4>w¹¹¹⁸, females had more overlapping genes responding to stress than males.

Separately, we analyzed 14 genes (that we used for PCR analysis), and found that some genes react to APS in a sex-specific manner. At the same time, genes whose expression strongly depends on the Gagr genotype are clearly identified: these are CG10232, CG10051, PPO1, and Nazo (Table 2).

Because males did not show a significant response to APS, we assessed the tissue specificity of the response only in females (Figure 5). As a result, we saw a systemic response covering a large number of tissues. In Tub-GAL4>w¹¹¹⁸ females, the largest number of genes that increased their expression in response to APS were expressed in the central nervous system, as well as in the digestive system. The largest number of genes that have decreased their expression were expressed in the gut, endocrine system, and reproductive system.

2.2.3. Transcription analysis of genes involved in immune pathways

Functional enrichment for genes that increased expression revealed many terms associated with stress response. Therefore, we separately assessed how APS affects the expression of genes involved in major stress signaling pathways.

We first assessed the JNK/MAPK and Jak-STAT stress cascades. It turned out that APS caused an increase in the expression of key transcription factors of the JNK cascade jra (Jun) and kay (Fos) in Tub-GAL4>w¹¹¹⁸ females (Figure 6). However, jra activation, although statistically significant, is below 50%. Similarly, statistically significant, but less than 50%, the transcription of the msn, Traf4 and Traf6 genes was activated. However, it was not found that any of the key kinases of the JNK cascade was regulated at the transcription level. Other JNK cascade components activated at the level of gene expression included Gadd45, whose function is associated with the regulation of the localization of JNK cascade proteins, and raw (encodes a membrane protein involved in dendrite patterning and the subcellular localization of JNK signaling components), as well as the puc gene, which encodes a serine/threonine protein phosphatase that mediates a feedback loop that regulates the Jun-N-terminal kinase pathway and forms negative feedback in the regulation of the JNK cascade [17]. The genes Pvf2 and Pvr, the ligand and its receptor activating the MAPK cascade, respectively, MAP kinase p38c also increased expression in response to APS; p38c MAP kinase encodes a protein involved in the stress and wound responses [18,19].

Thus, the JNK cascade was regulated at the level of expression of its components in response to APS. This regulation is associated exclusively with the activation of expression, extracellular ligands, their receptors and transcription factors, but not the main kinases.

In Tub-GAL4>w¹¹¹⁸ males, the transcription response of genes to APS was weaker. The p38c, ALiX, and jra genes were activated. The transcription response of genes to APS was poorly expressed in Tub-GAL4>UAS-Gagr males and females. Thus, knockdown of Gagr causes a disruption in the activation of the expression of genes involved in the JNK cascade.

We also analyzed the response to APS of the Jak-STAT signaling pathway genes (Figure 7). Exposure to APS leaded to activation of the expression of certain genes of the Jak-STAT signaling cascade in both Tub-GAL4>w¹¹¹⁸ females and males: upd2 (a ligand of the Jak-STAT cascade) and Socs36E (a negative regulator of the Jak-STAT cascade). However, transcription activation of the STAT92E transcription factor was observed only in Tub-GAL4>w¹¹¹⁸ females.

Thus, the Jak-STAT cascade is also regulated at the level of transcription of its components in response to APS. This regulation is associated with activation of the expression of cytokines, as well as a negative regulator of the cascade. Knockdown of Gagr causes a decrease in the activation of expression of the STAT92E gene in females, upd2, upd3 in males, but not the SocS36E gene (in both sexes) and upd2 in females.

We also examined the NFkB signaling pathways: Toll and Imd (Figure 8). Their role in innate immunity in Drosophila has been well studied [20], but little is known about their function in protecting the body when exposed to abiotic stress factors. In Tub-GAL4>w¹¹¹⁸ females and males, we did not find changes in the regulation of the expression of any of the Imd-signaling components in response to APS, except for a statistically significant, but less than 2 times, increase in the transcription of the genes ken in females and Dredd in males.

In turn, for the Toll signaling pathway, a change in the expression of several secreted factors (ligands, proteases, etc.) that positively regulate the activity of the Toll signaling pathway was detected: an increase in the expression of the gene GNBP2 coding the Gram-negative bacteria binding protein and the genes spz4, spz6 involved in Toll pathway-dependent AMPs production; a decrease in the expression of GNBP3 gene of the Gram-negative bacteria binding protein, the genes PGRP-SC1a, PGRP-SC2, PGRP-SD of peptidoglycan recognition proteins, the SPE gene coding a protease responsible for the cleavage of the Toll ligand. The Tub-GAL4>w¹¹¹⁸ males also responded to APS by decreasing in the expression of several genes.

Activation of the expression of intracellular Toll signaling components in response to APS was not detected in Tub-GAL4>UAS-Gagr males and females. Thus, knockdown of Gagr causes a disruption in the activation of the expression of genes involved in NfkB signaling pathways.

2.2.4. Genes whose expression is not induced by stress in the Gagr mutant

Next, we examined the number of genes that were not activated in Tub-GAL4>UAS-Gagr hybrids and that increased in expression in Tub-GAL4>w¹¹¹⁸ hybrids in response to stress. (Log₂FoldChange>1.5 for Tub-GAL4>w¹¹¹⁸ and Log₂FoldChange <0.5 for Tub-GAL4>UAS-Gagr). In females, there were 195 such genes. Functional category enrichment identified genes associated with development (Figure 9a). In males there were 115 such genes. Enrichment was shown for only one category – R-DME-975576: N-glycan antennae elongation in the medial/trans-Golgi.

Then, using the MetaScape tool, we carried out protein-protein interaction enrichment analysis for genes overexpressed in Tub-GAL4>w¹¹¹⁸ females (Figure 9b). A subset of representative terms was converted into a network layout. Each term is represented by a circle node, where its size is proportional to the number of input genes fall under that term, and its color represent its cluster identity (i.e., nodes of the same color belong to the same cluster). Terms with a similarity score > 0.3 are linked by an edge (the thickness of the edge represents the similarity score). The network is visualized with Cytoscape with “force-directed” layout and with edge bundled for clarity. Only physical interactions in STRING (physical score > 0.132) and BioGrid were used. The resultant network contains the subset of proteins that form physical interactions with at least one other member in the list. The circle size correlates with Log₁₀(P-value) from 10^-2 to 10^-10.

Enrichment of functional categories showed that genes activated in Tub-GAL4>w¹¹¹⁸ females and not activated in GAL4>UAS-Gagr females exhibit tissue specificity (Figure 9c), which is associated with expression in the crop, midgut, hindgut and central nervous system. A very similar pattern of tissue-specific expression is demonstrated by the Gagr gene (Figure 9d).

Since females demonstrated a more expressed pattern of activation of gene expression in response to APS, we next analyzed which genes were activated by stress in GAL4>w¹¹¹⁸ females and were not activated in the GAL4>UAS-Gagr females. Functional enrichment of a set of genes by molecular function showed that these genes form a group of 19 transcription factors (Fold enrichment 4.76; P value 8.51E-09, FDR 2.53E-05): run, ss, ase, sr, Antp, Sox21a, esg , grh, ham, Dfd, ich, nerfin-1, dmrt99B, grn, Kr-h1, acj6, rib, tap (Table 3). The biological functions of these gene products are primarily related to the development and functioning of the nervous system. Thus, in a mutant with knockdown of the Gagr gene, the disruption of many genes is apparently associated with impaired activation of a number of transcription factors. This means that their activation depends on the Gagr gene.

3. Discussion

We studied the transcriptomes of males and females with the Gagr gene knockdown in all tissues. In previous studies, we showed that this gene, derived from the gag gene of LTR retrotransposons, has a highly conserved structure across different Drosophila species, that indicate to a vital function of this gene. We also showed that this gene acquired a transmembrane domain during evolution and is most likely localized in the endoplasmic reticulum, as evidenced by protein-protein interactions of Gagr with the chaperone Pdi and the translation initiation factor eIF3j.

In the course of analyzing the transcriptome of the Gagr knockdown mutant, we found that in mutant flies, in comparison to the control strain of flies, the expression of immune response genes was increased (Figure 2 and Figure 3). These include the genes for antimicrobial peptides (AMPs) totC, AttA, AttB, AttC, AttD, DptB, LysB, and Dro, which suggest that the NFkB signaling pathway is not being properly regulated in the mutant. Furthermore, only the heat shock protein gene Hsp70Bb was upregulated.

It is known that activation of AMPs, including Drosocin, resulted in a significant extension of Drosophila lifespan [21]. AMP expression is regulated by members of the NFkB family, which include the transcription factors Dif, Relish and Dorsal. These transcription factors are activated upon infection by two major signaling cascades, the Toll and immune deficiency (Imd) pathways [22]. Additionally, subsets of AMPs can be directly activated by the transcription factors Drosophila Forkhead box O or Forkhead, depending on the metabolic status of the fly, demonstrating a cross regulation between metabolism and innate immunity [23]. In the midgut AMP expression is not regulated by Toll signaling but by the Imd and the Janus kinase-signal transducers and activators of transcription (Jak-STAT) pathways [24] and controlled by the negative transcription regulator caudal [25]. In our case, we did not find any significant changes in the expression of these transcription factors caudal, dfoxo, fkh. In addition, we did not find an increase in lifespan, but on the contrary, we recorded a reduced level of lifespan in flies with the Gagr knockdown (Figure 1a,b). This means that overexpression of AMPs by itself is not enough to increase lifespan and that in the mutant we studied, overexpression of AMPs is secondary and is a consequence of disruption of homeostasis processes.

However, we found an increased resistance of mutant females to heat stress (Figure 1c). The Turandot genes have shown promise as a mechanism that may contribute to inducible thermal tolerance and explain the asymmetry between expression level of Heat shock protein, and improved thermal tolerance. In the study [26], it was found that Drosophila with an overexpression of TotC and hsp70 was functionally linked to increased heat tolerance; at the same time expression of totA did not have a significant effect on heat stress tolerance [27]. Thus, the increased resistance of the Gagr gene mutant to heat stress may be associated with increased expression of the AMP genes and the Hsp70Bb gene.

We found melanotic nodules in mutant females (Figure 1c). Apparently, this may indicate the induction of an autoimmune reaction. It is known that the humoral response leads to the synthesis of dedicated AMPs and contributes to hemolymph coagulation and melanization. Crystal blood cells are involved in melanization. Crystal cells express the enzyme phenoloxidase, responsible for the initiation of the melanogenesis cascade. The activation of Pro-phenoloxidase is partially controlled by the serine protease inhibitor serpin 27A (Spn27A). Spn27A mutant larvae shows a melanotic phenotype and excessive melanization in response to immune challenge [28]. This phenotype is linked to the activation of the Toll pathway [29,30]. The Toll pathway controls the nuclear targeting of the Drosophila NF-κB proteins Dorsal and Dif. Constitutive activation of the pathway in Toll gain-of-function or cactus loss-of-function mutants leads to overproliferation of hemocytes, in particular lamellocytes, resulting in the formation of melanotic nodules [31,32].

The other signaling pathways also can activate melanozation process. After immune challenge of larvae with constitutive expression of the PGRP-LE, functioning upstream of the Imd pathway, melanotic masses were observed in the cuticle and hemolymph [33]. The activation of other pathways like the Ras/MAPK in hemocytes by the expression of transgenes leads to hemocyte proliferation and formation of melanotic masses [34,35]. Constitutive activation of Jak-STAT signaling resulting from a dominant Jak mutation, hop^Tum-l, induces the TotA gene upregulation, and an overproliferation of circulating plasmatocytes and differentiation of lamellocytes, leading to the formation of melanotic masses in larvae and adult flies [36]. Also, it was demonstrated that tumor(1)Suzukit (tuSz¹) mutant displayed a temperature-sensitive self-encapsulation phenotype directed at its own posterior fat body tissue [37,38]. Recently, it has been shown, that the tuSz¹ mutant phenotype may be the result of a gain-of-function mutation in the hop gene, but also a loss of function mutation in the GCS1 gene, which disrupts the protein N-glycosylation pathway in the posterior fat body [39]. The obtained results demonstrate that N-glycosylated extracellular matrix proteins serve as self-associated molecular patterns (SAMPs) and that activated innate immune cells attack tissues that lack these SAMPs. The self-tolerance mechanism can also act as an initiator of immunity itself if failure to recognize a self-signal on the pathogen surface triggers an immune response, in what is known as “missing-self recognition” [40].

It should be noted that the transcription of the GCS1 gene in both GAL4>UAS-Gagr females and males, is not significantly changed. Thus, it is impossible to associate melanization in females with its function. However, under stress conditions, GCS1 expression is significantly reduced in Gagr females (Fog₂FoldChange=-0.59, Padj=0.0001).

We observed melanization only in the Gagr mutant females. Of the three genes of the Ppo family, only the Ppo1 gene showed statistically significant changes in the transcription level. Moreover, the expression of the Ppo1 gene in both mutant and control flies is lower in females than in males. At the same time, in flies with the Gagr gene knockdown, the expression of the Ppo1 gene is lower relative to the control flies, and decreases upon induction of stress (Fog₂FoldChange=-0.84, Padj=0.02). The apparent contradiction is can be explained by the fact that we measure the level of gene transcription throughout the body, but not in crystal cells. According to FlyBase, Ppo1 expression is observed at high levels in muscle cells and carcass, and overall expression is higher in males than in females. Our results suggest that the mutant undergoes a signal switch that results in decreased expression of this gene throughout the body, and likely site-specific activation in crystal cells.

We found that transcriptomic response to stress in males is less significant than in females (Figure 4 and Figure 5, Table 2). Sexual dimorphism in the immune response is well known and characterized [41]. It should be noted that the transcription of the Gagr gene is characterized by sexual dimorphism: it is expressed at a level approximately in 2 times higher in males than in females; and its expression is not induced by APS [8]. Apparently, the low level of gene activation in response to stress in the mutant males is explained by a higher level of expression of immune response genes under standard conditions.

Since our previous studies have shown that transcription of the Gagr gene is most noticeably induced only in females by a strong oxidant, ammonium persulfate. Therefore, in experiments to study the transcriptomic response to stress, we used this agent. APS exerts its effect mainly outside the cell and acts on membrane proteins on cell surface, which usually gives rise to a decrease in cell viability and an increase in cell apoptosis [42]. APS produced a high level of oxidative stress in lysosome and induced epithelial-mesenchymal transition (EMT) via lysosomal oxidative stress [43]. Thus, APS causes a powerful stress response. However, studies of the transcriptomic response to APS in Drosophila have not yet been conducted.

We found that APS causes significant changes in the transcriptome of females, causing the activation of genes associated with protective stress responses (activation of the immune response, inflammation, chitin metabolism), and the suppression of genes involved in the metabolism of fats, proteins and carbohydrates and genes involved in oogenesis (Figure 3). Components of stress-associated signaling pathways JNK, Jak-STAT, Toll are regulated at the transcription level in response to APS. This regulation is mainly associated with transcription activation. The transcriptomic response to APS occurs in all organs; activation of expression in the digestive system is potentially regulated by the transcription factors kayak (JNK cascade) and Stat92E (Jak-STAT signaling pathway), in the fat body and carcass – by Stat92E, and in the head – by the transcription factor of the Toll pathway, Dif (Figure 5). It is noteworthy that the tissue specificity of the response to APS coincides with the tissue specificity of Gagr gene transcription (Figure 9).

Knockdown of the Gagr gene leads to disruption of the normal activation of stress-associated signaling cascades - JNK, Jak-STAT and Toll (Figure 6, Figure 7 and Figure 8). The transcriptomic response to APS is weakly expressed. Transcription of genes that are activated in the control flies and not activated in the Gagr gene mutant flies is associated with the activity of genes in the digestive system and central nervous system (Figure 8, Table 3). It is noteworthy that the tissue specificity of the Tub-GAL4>w¹¹¹⁸ response to APS coincides with the tissue specificity of the transcription of the Tub-GAL4>UAS-Gagr gene and the Gagr gene (Figure 9). Overexpressed genes were enriched in one category—transcription factor genes, whose function is associated with the regulation of development and functioning in the central nervous system.

In males, genes associated with N-glycosylation were significantly activated in the control flies and not activated in the mutant. Notably, N-glycosylation is associated with immune responses. It is known that many immune proteins and antibodies are glycosylated. Pathogen's glycoproteins play vital roles during the infection cycle and their expression of specific oligosaccharides via the N-glycosylation pathway to evade detection by the host immune system [44]. N-glycans are important for the processes that precede or follow the actual sorting event, such as protein folding, quality control, endoplasmic reticulum (ER)-associated degradation, ER-to-Golgi trafficking, and retention of glycoproteins in the apical membrane [45].

The Gagr gene is closely integrated into the regulatory network of signaling cascades: its transcription depends on signals from the JNK and Jak-STAT pathways [8]. The assumption of this regulation is in good agreement with the results of our experiments and other studies in which the activation of Gagr expression is observed in response to significant stress (viral infection, oxidative stress caused by peroxo compounds). The JNK signaling pathway has many functions, regulating a diversity of processes from cell movement during embryogenesis to the stress response of cells after environmental insults. In both Drosophila and higher organisms, JNK takes part in different processes, including apoptosis, proliferation, differentiation, cell migration, tumorigenesis, cell competition, and processes of cell regeneration [46,47]. The kayak protein (a part of the AP-1 transcription factor) is a developmentally regulated transcription factor that may play a role in the function or determination of a particular subset of cells in a developing embryo [48]. In wounded tissues, JNK is activated in the damaged cells to ensure their apoptotic death and in the surviving cells to promote their cellular reprogramming and proliferation [49].

It is known that JNK and Jak/STAT activation in imago promote the proliferation of stem cells (SCs) in response to oxidative or ER stress and infection [50]. In addition, the JNK pathway regulates upd3 (effector of Jak/STAT pathway) expression, which is necessary for optimal renewal of the intestinal epithelium and survival following septic injury [48]. JNK also becomes widely activated in the intestinal epithelium of aging flies, inducing excessive proliferation of ISCs [51]. In addition, autophagy plays a role in Drosophila ISCs to maintain proliferation and preserve the stem cell pool.

It can be concluded that that the expression of the Gagr gene is observed in imago tissues with a high potential for stress-induced proliferative activity. Thus, Gagr likely participates in the control of morphogenesis at the embryonic stage of development, and in adults, in post-stress tissue regeneration. Based on the proven localization of the Gagr protein as well as some of its partners in membrane, we can conclude that it is inserted into the ER membrane and binds to a component of the translation system and signaling proteins. Knockdown of the Gagr gene likely leads to global changes in gene expression, confirming its important role in maintaining homeostasis. The next step in our investigation is tissue specific transcriptome analysis.

4. Materials and Methods

4.1. Drosophila melanogaster strains and conditions

The following strains of D. melanogaster were used in the work: w¹¹¹⁸; driver Tub-GAL4 (y[1] w[1118]; P{w[+mC]=tubP-GAL4}LL7 P{ry[+t7.2] =neoFRT}82B/TM6B, Tb[1]), the expression of GAL4 in which is controlled by the tubulin gene promoter; VDRC 107457 KK, carrying a transgenic construct for expression of dsRNA for the Gagr gene RNAi under UAS control (hereinafter referred to as UAS-Gagr). Strains were obtained from the Vienna Drosophila RNAi Center. Fly stocks were maintained in a standard nutrient agar medium at 25°C. To induce interference, females of the UAS-Gagr strain were crossed with males of the Tub-GAL4 driver strain. As a control, females of the w¹¹¹⁸ strain were crossed with Tub-GAL4 males. Adult flies 6–8 days old were used in all experiments.

4.2. Physiological tests

To analyze lifespan, flies at the age of 1 day were selected, separated by sex, and put into separate test tubes of 20-30 flies. The number of individuals in the test tube was checked every 2-3 days, the food was replaced every 5 days. Lifespan was measured at 27°C in medium supplemented with 0.1 M APS and standard medium. To measure the mobility of adults, the climbing test was used [52]: 30 adults were placed in an empty long tube. Flies were dropped to the bottom of the test tube by mechanical tapping. Next, the maximum vertical rate of rise to a height of 17 cm was assessed. Two repetitions were performed. Heat shock survival study was conducted in a water bath at 38°C. The flies were placed in a water bath for 2 hours, after which they were left in a thermostat at 27˚C for 24 hours. After 24 hours, the number of dead flies was counted. Two repeats were made: for the Tub-GAL4>UAS-Gagr hybrids there were 24 flies in the first repeat, 29 flies in the second, and for Tub-GAL4>w¹¹¹⁸ there were 23 flies in the first repeat, and 27 flies in the second.

4.3. RNA extraction and RT-PCR

Total RNA was isolated from Tub-GAL4>UAS-Gagr and Tub-GAL4>w¹¹¹⁸ flies after 24-hour exposure to 0.1 M APS in females and males aged 5-7 days. Flies reared under standard conditions served as controls. RNA was isolated from pools (five females, seven males) in 5–7 biological replicates using the ExtractRNA reagent (Evrogen, Moscow, Russia), according to the manufacturer's protocol, then it was treated with DNase I (Thermo Fisher Scientific, USA). Reverse transcription was carried out using an MMLV-RT kit (Evrogen, Russia), according to the manufacturer’s protocol, with random primers (Evrogen, Russia). For quantitative PCR with the obtained cDNA, a Taq polymerase-based reaction mixture with intercalating dye SYBR Green I (Evrogen, Russia) was used in accordance with the manufacturer’s protocol. The reaction was performed using a Mini Opticon Real-Time PCR System (Bio-Rad Laboratories, Hercules, CA, USA). The relative expression of the genes Gagr, TotA, TotC, AttB, CecA2, Socs16D, Spn88Eb, CG33346, Nazo, Ppo1, Spn28Dc, CG1304, CG10232, Ser6, CG10051, normalized to the expression of two reference genes, Tub84D and EloB, was analyzed. Amplification was performed with primers shown in Table S3. To present the expression results, a histogram was constructed in the GraphPad Prism 9 program (https://www.graphpad.com/). Statistical significance was assessed using the nonparametric Mann-Whitney test in GraphPad Prism 9.

4.4. RNA-sequencing and data processing

Total RNA was isolated from UAS-Tub-GAL4>UAS-Gagr and Tub-GAL4>w¹¹¹⁸ flies after 24-hour exposure to 0.1 M APS in females and males aged 5-7 days. Flies maintained under standard conditions without APS served as control. RNA was isolated from pools (five females, seven males) in 5–7 biological replicates using the ExtractRNA reagent (Evrogen, Moscow, Russia), according to the manufacturer's protocol, then it was treated with DNase I (Thermo Fisher Scientific, USA). RNA concentration and integrity were evaluated by a fluorimetric assay (Qubit 4, ThermoScientific, USA) and capillary electrophoresis (TapeStation, Agilent, Germany), respectively. All samples were prepared in one experiment (3 repeats for each sample). Strand-specific libraries were prepared by the NEBNext Ultra II Directional RNA Library Preparation kit (NEB, USA) and sequenced (100 nucleotides, single end) with median depth 25 million reads per sample by NovaSeq 6000 (Illumina, USA). Low quality reads and adapter sequences were deleted (Timmomatic tool, v0.36), then reads were aligned to the BDGP6 primary genome assembly. Uniquely aligned reads were counted for known exons of each gene using R subread package (R environment). For the reference genes ppl, Tbp, Gapdh1, tub, RPL40, SdhA, the expression deviation for each gene in the sample was assessed, normalized to the expression value of the corresponding gene in the control samples without exposure to APS. Deviations for different genes averaged +/-0.3. Differential expression analysis performed by DESeq2 package (R) (https://bioconductor.org/packages/release/bioc/html/DESeq2.html). Differentially expressed genes (DEGs) were defined as protein coding genes with Padj<0.05, |Fold Change|≥1.5. To analyze overlapping gene samples, Venn diagrams were constructed (https://bioinformatics.psb.ugent.be/webtools/Venn/). Metascape analysis tools were used to identify functional enrichment categories of DEGs (http://metascape.org; [53]). GeneOntology Resource (https://geneontology.org/) was used to search for molecular function enrichment. In order to determine which tissues the transcriptome response detected at the level of the whole organism may be associated with, we used data from the FlyAtlas2project [54]. We assessed correlations of tissue-specific expression values for differentially expressed genes and built a gene co-expression network for correlations of 0.8 or more using Cytoscape using the graph-oriented clustering method MCODE (https://cytoscape.org/).