text:
Construction and Analysis of the Complete Genome Sequence of Leprosy Agent Mycobacterium lepromatosis

ABSTRACT Leprosy is caused by Mycobacterium leprae and Mycobacterium lepromatosis. We report construction and analyses of the complete genome sequence of M. lepromatosis FJ924. The genome contained 3,271,694 nucleotides to encode 1,789 functional genes and 1,564 pseudogenes. It shared 1,420 genes and 885 pseudogenes (71.4%) with M. leprae but differed in 1,281 genes and pseudogenes (28.6%). In phylogeny, the leprosy bacilli started from a most recent common ancestor (MRCA) that diverged ~30 million years ago (Mya) from environmental organism Mycobacterium haemophilum. The MRCA then underwent reductive evolution with pseudogenization, gene loss, and chromosomal rearrangements. Analysis of the shared pseudogenes estimated the pseudogenization event ~14 Mya, shortly before species bifurcation. Afterwards, genomic changes occurred to lesser extent in each species. Like M. leprae, four major types of highly repetitive sequences were detected in M. lepromatosis, contributing to chromosomal rearrangements within and after MRCA. Variations in genes and copy numbers were noted, such as three copies of the gene encoding bifunctional diguanylate cyclase/phosphodiesterase in M. lepromatosis, but single copy in M. leprae; 6 genes encoding the TetR family transcriptional regulators in M. lepromatosis, but 11 such genes in M. leprae; presence of hemW gene in M. lepromatosis, but absence in M. leprae; and others. These variations likely aid unique pathogenesis, such as diffuse lepromatous leprosy associated with M. lepromatosis, while the shared genomic features should explain the common pathogenesis of dermatitis and neuritis in leprosy. Together, these findings and the genomic data of M. lepromatosis may facilitate future research and care for leprosy. IMPORTANCE Leprosy is a dreaded infection that still affects millions of people worldwide. Mycobacterium lepromatosis is a recently recognized cause in addition to the well-known Mycobacterium leprae. M. lepromatosis is likely specific for diffuse lepromatous leprosy, a severe form of the infection and endemic in Mexico. This study constructed and annotated the complete genome sequence of M. lepromatosis FJ924 and performed comparative genomic analyses with related mycobacteria. The results afford new and refined insights into the genome size, gene repertoire, pseudogenes, phylogenomic relationship, genome organization and plasticity, process and timing of reductive evolution, and genetic and proteomic basis for pathogenesis. The availability of the complete M. lepromatosis genome may prove to be useful for future research and care for the infection.

RESULTS AND DISCUSSION
Construction of the complete genome. The construction of M. lepromatosis genome encountered several technical hurdles. Due to the lack of axenic cultivation, the quantity and quality of genomic DNA became the first rate-limiting step. The genomic DNA was extracted from a dried, stained, and archived smear slide, with only ;18% of sequenced reads (12.3 of 69 million reads) belonging to M. lepromatosis. Contaminant reads included human host DNA and nonspecific bacteria DNA. Some bacterial DNA reads, from Thermus scotoductus particularly, were likely introduced and/or amplified during the library preparation step that used polymerase chain reactions (PCR) to boost target quantity. Other nonspecific bacterial DNA, such as from Propionibacterium acnes, Micrococcus luteus, and others, were likely from the stained smear that used regular reagents for acid-fast stain, tap water, and microscopy immersion oil. To overcome contamination, initial steps used matches to the nearest M. leprae genome to capture specific reads, de novo assembly, systematic tagging and removal of contaminant contigs, and alignment to M. leprae for draft construction. Later more iterative cycles of mapping, de novo assembly, and gap filling were used to retrieve unique sequences, close gaps and refine drafts. PCR and Sanger sequencing of amplicons were used to close final gaps, verify assembly of several long contigs, and settle copy numbers of short repetitive sequences. The genome coverage was 525Â.
A phylogenomic tree was constructed based on the proteomes (316,716 alignment sites) of well-known mycobacteria and the newly named and sequenced Mycobacterium uberis, a yet-to-be cultivated organism that causes nodular thelitis and tuberculoid scrotitis in cows and goats (22). As shown in Fig. 1, M. leprae and M. lepromatosis stem from a most recent common ancestor (MRCA) that forms a clade with M. haemophilum and M. uberis. Shortly after the divergence of M. haemophilum, the lineages of M. uberis and the MRCA diverged again and both underwent genome downsizing, pseudogenization, and decrease in GC content (Table 1).
Comparison with the draft genome of M. lepromatosis Mx1-22A. The FJ924 genome was compared with the draft genome of M. lepromatosis Mx1-22A that consisted of 3,206,741 bp from 126 contigs (10). This detected 66 kb extra nucleotides from 56 segments with lengths .300 bp, with the largest segment covering 8,518 bp (Table  S1A). Annotation of these segments revealed 27 CDS (21 kb, 1.5% of all CDS) and 65 pseudogenes (44 kb, 4.2% of all pseudogenes). Most of these genes and pseudogenes were absent in Mx1-22A, while a few of them showed partial sequence of the  3  3  3  3  3  2  3  tRNA genes  46  45  45  45  45  44  45  tmRNA gene  1  1  1  1  -1  1  Other noncoding RNA genes  13  3  3  2  3  3 (Table S1B and S1C). Shorter Mx1-22A was likely due to loss of some M. lepromatosis unique sequences during array capture to the M. leprae genome, to the assembler that discards/collapses repetitive regions, and to low sequencing coverage of the genome (55Â) (10). Single nucleotide polymorphisms (SNPs) between Mx1-22A and FJ924 were identified with MAUVE (23,24). After exclusion of ambiguities and inaccuracies near the edges of Mx1-22A contigs, a final set of 11 SNPs were obtained, with six in genes (five in CDS and one in the 16S rRNA gene) and five in intergenic regions (Table S2). Four SNPs detected in coding genes produced one amino acid difference between the proteins encoded by the two strains. The fifth SNP in a CDS encoded a PPE protein of 410 amino acids in FJ924 (MLPF_1455), but it caused a stop codon in Mx1-22A, hence a pseudogene (MLPM_1054). The limited number of SNPs suggests that the divergence between Mx1-22A and FJ924, strains likely from the Monterrey region in north central Mexico and the Sonora region in north western Mexico, respectively (;1000 km away) (4,10), was very recent, probably a few hundred years ago. Similarly, another M. lepromatosis strain PL-02, also from Mexico, diverged from Mx1-22A ;186 years ago (21).
Among the extra CDS detected, the complete FJ924 genome contained two tandemly arranged asparagine permease genes (MLPF_1849 and 1850), like M. leprae, and three identical copies of the gene encoding bifunctional diguanylate cyclase/phosphodiesterase, unlike M. leprae, M.uberis and M. tuberculosis with a single copy (further analysis later). The Mx1-22A contained single asparagine permease gene and single diguanylate cyclase/phosphodiesterase gene (10); such omissions were likely due to relaxed genome assembly and/or lack of copy count instead of true absence. Assembly errors in three Mx1-22A contigs were also seen, such as contig-111, a 14,151-bp fragment, that showed a minus orientation of the first 8,706 bp.
Repetitive sequences. The M. leprae genome is characterized by the presence of numerous repetitive sequences that likely play a role in genome plasticity (25). Repetitive sequences that appeared twice or more in the M. lepromatosis FJ924 genome were analyzed (Table S3). In total, they accounted for 3.5% of the genome and were scattered in 125 genome segments. Some repeats dispersed similarly to the repetitive elements in the M. leprae genome (25). Four main families of repetitive elements were identified: types A, B, C and D, with 38, 8, 8 and 6 copies, respectively ( Fig. 2 and Table S3A). These repeats were named as RLPM, LPMREP, REPLPM, and LPMRPT to correspond to RLEP, LEPREP, REPLEP, and LEPRPT, respectively, in the M. leprae genome.
Among the 38 copies of RLPM of various extension lengths, a core sequence of 626 nucleotides was identified in 18 copies (see consensus in Table S3B). The complete or partial The Complete Genome of M. lepromatosis Microbiology Spectrum core sequences were highly conserved in the 38 copies (99.2-100% identity). However, compared with all 37 RLEP copies in M. leprae, RLPM showed low sequence identity, in 25 copies at best 75% of 126 nucleotides with several gap sites. Nonetheless, in view of the nearly identical copy numbers, similar average length (;700 bp), and some homologous genome locations, RLEP and RLPM probably evolved from the same ancestral repeat after divergence from MRCA. Type B repeats LPMREP resembled the LEPREP repeats in M. leprae in BLAST matches (79-83%) and copy numbers. Both showed five complete copies and three fragment copies in their respective genomes. The origin of these repeats seems to be a pseudogene of a putative group II intron maturase-related protein, which persisted from the MRCA. Based on the alignment of the complete copies of LPMREP, a core sequence of 2,453 nucleotides  (Table S3B). Recently, these repeats were used in a M. lepromatosis molecular diagnostic assay, but under the abbreviation of 'RLPM' (16). Type C repeats REPLPM in M. lepromatosis were equivalent to the REPLEP repeats in M. leprae, but with gapped alignments and an overall 74% match. The copy numbers also varied: 8 copies for REPLPM but 15 copies for REPLEP (25). The REPLPM copies displayed different lengths but a consensus core of 1,123 nucleotides could be inferred (based on 50% coverage). This core sequence (Table S3B) was completely included in REPLPM_07.
Type D repeats LPMRPT were equivalent to LEPRPT repeats in M. leprae with 87% similarity covering almost the complete sequence. Only four of LPMRPT repeats were large (;1,200 bp) (see LPMRPT_01 in Table S3B), while the two others were short fragments.
Functionally, the four repetitive elements likely also played a role in genome plasticity and organization (see below). Other repetitive sequences were also annotated (Table S3A). They likely arose from segmental duplications or gene families.
Comparative genomics: synteny and proteome. Synteny among the genomes of M. leprae, M. lepromatosis, M. haemophilum and M. tuberculosis was compared with the program MAUVE (23) (Fig. 3). Many differences were observed between the two leprosy bacilli, by one way, and between M. haemophilum and M. tuberculosis, by the other, indicating events of chromosomal rearrangement.
To infer the minimal number of chromosomal rearrangements among them, a pairwise inversion distance matrix was constructed with the program GRIMM (26), which predicted 18 inversion events between M. leprae and M. lepromatosis and 11 between M. haemophilum and M. tuberculosis (Fig. 3). In addition, the MRCA of leprosy bacilli incurred ;50 chromosomal rearrangements after divergence from M. haemophilum. Comparisons of the consecutive CDS and pseudogenes in M. leprae and M. lepromatosis further rendered 25 syntenic blocks as the likely outcome of major rearrangement events (Table S4). These blocks and the 24 synteny breaks were analyzed for genesis by comparing them with ancestral M. haemophilum and M. tuberculosis (Table S5) (Table S5).
At the 24 synteny breaks, 17 repetitive elements were noted to be involved: 9 RLPM, 3 REPLPM, 2 LPMREP and 3 other repeats ( Fig. 2 and Table S5). These observations reveal the abrupt effect of repeats through recombination to shape the genome of M. lepromatosis. These chromosomal rearrangements may also produce gene duplications and deletions and affect gene expression. Similar effects were also previously observed for the genome of M. leprae (25).
Orthologous clusters of proteins were estimated with the program OrthoFinder v2.3.3 (27) for proteomes of the four mycobacteria (Fig. 4). A core proteome of 1,332 The relative ages of pseudogenes (Rpa) are shown in Fig. 5. Values for unique pseudogenes were small, with mean and SD of 0.138 6 0.139 and 0.216 6 0.190 for M. lepromatosis and M. leprae, respectively, suggesting more recent occurrences after divergence. Values for shared pseudogenes were larger, 0.462 6 0.238 for M. lepromatosis and 0.467 6 0.254 for M. leprae along with a wider and overlapping distribution. An xy plot of the shared pseudogenes showed a positive correlation (R 2 = 0.704) (data not shown), indicating a common pseudogenization process prior to divergence. Together, the analysis of Rpa suggests that most of them occurred approximately during the same period as an event of massive pseudogenization.
A time tree analysis using the RelTime method (28) estimated the time of divergence between M. haemophilum and the MRCA of the leprosy bacilli to be 29.52 Mya (Fig. 6). When this value was set to be Rpa of 1, Mya values for the shared pseudogenes  (10), but much later than the divergence from the related species M. uberis at 28 Mya (Fig. 6).
Pathogenetic insights from shared genomic features. The leprosy bacilli have a moderate pathogenicity, between rarely pathogenic M. haemophilum and highly pathogenic M. tuberculosis (11). The lean genomes and shared genomic features underscore the clinical and pathological features of leprosy, i.e., long incubation period, insidious onset, dermal and neuronal invasion of bacillus-laden macrophages, and chronicity, contrasting to mainly pulmonary infection in tuberculosis. Table 2 shows the numbers of shared and unique CDS and pseudogenes for M. leprae and M. lepromatosis. The two bacilli shared 1,420 CDS and 885 pseudogenes (71.4% of genome) but differed in 1,281 CDS and pseudogenes (28.6%), although some of them may be associated with the different annotation protocols. Some shared features were examined to highlight common pathogenesis in leprosy.
Among all known mycobacteria (;160 species), neuronal invasion is unique to the leprosy bacilli. In M. leprae, adherence to Schwann cells seems to require the presence of the DNA-binding protein HU and the cell wall antigen phenolic glycolipid 1 (PGL-I) (29). PGL-I also subverts host immune cells to shield the bacillus from clearance (30). The synthesis of PGL-I requires 6 enzymes/genes (31). In M. lepromatosis FJ924, the PGL-I genes were identified as MLPF_0154, 0155, and 0156, and MLPF_2688, 2689, and 2690 in two clusters without related pseudogenes. The bacillus also contained the gene encoding DNA-binding protein HU (MLPF_0968) that is also known as mycobacterial DNA binding protein 1 or laminin-binding protein (32). This protein consists of an N-terminal HU-like region and a C-terminal IDR (intrinsically disordered region). The M. lepromatosis protein contained 195 residues, slightly smaller than the proteins of M. leprae, M. haemophilum and M. tuberculosis (200, 213 and 214 residues, respectively). The length difference was close to the C terminus at the IDR, essential region for such function as genome compaction or suppression of DNA synthesis (33). Several ESX secretion systems (type VII secretion systems or early secretory antigenic target [ESAT6] secretion) are important for virulence of mycobacterial species (34,35). Among them, ESX-1 system is the main determinant of virulence in M. tuberculosis. The genes encoding the proteins of the ESX-1 system are located in two or three clusters in the genomes of M. tuberculosis, M. haemophilum and M. uberis (22). In M. lepromatosis Mx1-22A and M. leprae, around a quarter of the genes are nonfunctional pseudogenes (10). In M. lepromatosis FJ924, the genes and pseudogenes of the ESX-1 system were located at two clusters in the genome, one from pseudogene MLPF_0041 (espE) to gene MLPF_0062 (mycP) (one of the three mycP genes of M. lepromatosis) and the other comprising genes MLPF_601, 602, and 603 (espA, espC and espD). Retention of some ESX-1 genes in the leprosy bacilli suggests some functionality of the secretion system.
The ESX-5 system plays a pivotal role in the secretion of proteins rich in Pro-Glu (PE proteins) and Pro-Pro-Glu (PPE proteins) (35). M. lepromatosis contained mainly the core genes of the system, MLPF_0876 to 0880 (eccA, eccE, mycP, eccD and espG) and MLPF_0889 and 0890 (eccC and eccB). A pseudogene of the esat-6 like protein (MLPF_0881) was detected between these two gene clusters. In the orthologous M. leprae region, neither an esat-6 gene nor a pseudogene was detected. Genes encoding the PE/PPE proteins were examined in the genomes of leprosy bacilli for comparison with those in M. haemophilum and M. tuberculosis. OrthoFinder clustered these proteins/genes in 10 orthogroups (Table S6). While the M. haemophilum and M. tuberculosis genomes each contained 46 genes spread in the clusters, M. lepromatosis and M. leprae each contained only 13 genes but 34 and 32 pseudogenes, respectively (Table S7).  PE/PPE proteins are noted to be antigenic determinants in M. tuberculosis to play a role in immuno-pathogenesis (36). To the less virulent leprosy bacilli, these alarming molecules would be detrimental to the survival. Rather, in M. leprae, preferential inactivation of PE/PPE genes is noted whereas PGL-I and its genes, beneficial to the bacillus, are retained, which likely aid immune evasion and survival of the bacillus during ;20 million years of reductive evolution (11). In M. lepromatosis, these were also true, with statistical significance ( Table 3). As a benchmark for gene retention, all 52 genes but one that encode ribosomal proteins, essential for ribosome structure, protein synthesis and bacterial survival, remained ( Table 3). The same was true for M. leprae. The similar numbers of genes and/or pseudogenes for PE/PPE, PGL-I, and ribosomal proteins suggest their origins in the MRCA and persistence.
Dermal tropism is innate to the clade of M. haemophilum, lineages of leprosy bacilli, and M. uberis. Despite unknown environmental niches, M. haemophilum grows optimally at 30-32°C in culture and requires heme supplement (12,37). The bacillus causes rare infections, mainly dermal and/or disseminated ones, in immunocompromised patients, such as long-term steroid users. In patients with leprosy, the most common sites of infection involve earlobes, extremities, nasal mucosa, and scrotum where the temperatures tend to be below 37°C to favor bacterial proliferation.
Finally, the large number of shared genes and pseudogenes would require a systematic approach to analyze their functional clusters and/or lack thereof. In view of the cultivation difficulty that has impeded the research for leprosy bacilli, future research endeavors in M. tuberculosis and M. haemophilum should help in this regard.
Pathogenetic insights from unique genome features. The two leprosy bacilli differed in CDS and pseudogenes, including 369 unique CDS in M. lepromatosis and 184 unique CDS in M. leprae ( Table 2). Examination of these unique CDS and pseudogenes should offer insights into some varying pathological features in patients with each infection, such as more skin nodules associated with M. leprae but more invasion into subcutis, vessels and internal organs associated with M. lepromatosis.
The finding of three copies of the diguanylate cyclase/phosphodiesterase gene in M. lepromatosis (MLPF_1053, 2468 and 2944) contrasted the single copy status in M. leprae (MLBr01750), M. uberis (WP_116539319.1) and M. tuberculosis (Rv1354c) but resembled those in M. haemophilum. A phylogenetic reconstruction of the evolution of these proteins in several mycobacterial species (Fig. 7) revealed that the encoding gene was duplicated in the early ancestor of these species and Mycobacterium avium (taxonomy in Fig. 1). After this event, several gene gains and losses took place. M. lepromatosis, M. leprae and M. uberis retained the same copy and lost the other, in opposition to M. tuberculosis. Gene MLPF_1053 was located at a large syntenic block of . 200 kb (block 11, Table S4), suggesting ancestral origin. Gene MLPF_2468 arose through a 5,020-bp duplication involving MLPF_1053 and 1052 (a pseudogene derived from the transposase gene of an insertion sequence element) and the last part of pseudogene MLPF_1051. Gene MLPF_2944 was a duplicate of a 4,700-bp segment involving MLPF_2468 and the pseudogenes MLPF_2467 and 2466. In M. leprae, the amplified copies as well as associated pseudogenes were absent (Table S4). Surprisingly, in a recent M. leprae genome, MRHRU-235-G from India (CP029543.1, direct submission), a tandem gene duplication was also noted. The diguanylate cyclase/phosphodiesterase is a well-studied enzyme in diverse bacteria, harboring three domains: GAF as sensor domain, GGDEF as diguanylate cyclase, and EAL as c-di-GMP specific phosphodiesterase (38). The latter are two antagonistic enzymes, hence bifunctional, involved in the synthesis and hydrolysis of the second messenger cyclic di-GMP. The presence of three copies of its gene hints more c-di-GMP metabolism in M. lepromatosis for survival and/or pathogenicity in view that deletion of this gene in M. tuberculosis affects its dormancy and pathogenicity (39,40). In M. leprae, studies have also suggested a potential role of diguanylate cyclase in the signaling response for intracellular survival (41). In other pathogens, a strategy has been designed to intercept c-di-GMP signaling pathways by directly targeting this second messenger as a way of controlling pathogen survival, growth or biofilm formation (42).
The gene hemW, previously annotated as hemN (10), was found in M. lepromatosis (MLPF_2234) but not in M. leprae. This gene likely encodes a heme chaperone that catalyzes the insertion of heme into hemoproteins for respiration (43,44). The gene and adjacent MLPF_2235 (a pseudogene) stood between syntenic blocks 19 and 20 (Table  S4); in M. leprae, however, the syntenic blocks rearranged afar after divergence, during which hemW and adjacent pseudogene (together ;2 kb) were likely lost. In M. haemophilum and M. tuberculosis, this gene was B586_09185 and Rv2388c, respectively, and their corresponding amino acids shared 89% and 81% identities with M. lepromatosis. This finding warrants further research into the function of hemW.
The two leprosy bacilli showed 6 genes in common to encode TetR family transcriptional regulators (ethR and others), and there were 5 additional ones in M. leprae but only pseudogenes in M. lepromatosis. These regulators are involved in diverse functions, such as regulation of c-di-GMP, ethionamide resistance, lipogenesis, etc. (45,46). The difference in these regulators, along with many other unique genes, suggest considerable variations between the leprosy bacilli in metabolism, stress response, and interaction with the host to influence pathogenesis.
Conclusions. The history of the leprosy bacterial lineage started around 30 Mya when a mycobacterial ancestor diverged from the lineage of M. haemophilum. At that time, this ancestor would have a similar lifestyle to the present bacterium M. haemophilum, perhaps being able to produce opportunistic skin infections (12). Soon after this point, two lineages evolved, one leading to M. uberis, a pathogen infecting cows and goats (22) and the other leading to the MRCA of M. leprae and M. lepromatosis. Based on the phylogenies and pseudogene comparisons, the lineages of M. uberis and the MRCA likely underwent independent processes of reductive genome evolution (22). While the two leprosy lineages diverged around 13.9 Mya (10), our estimation of the ages of pseudogenes suggest that most of them formed shortly before the divergence. Afterwards, additional events of pseudogenization took place independently in both lineages, but at a much lesser scale. We proposed previously that the reductive evolution of leprosy bacilli began in the MRCA as a parasite within the host, family Hominidae (great apes), ;18 to 20 Mya (11). In another word, the specific taming-adapting process has persisted all along until modern humans today. Our current genomic study refines the bacterial reductive evolution to ;14 Mya, which brings the host more specifically within subfamily Homininae (humans, gorillas, and chimpanzees). The genome reduction of M. uberis that also signifies a host specific process lends further support. Likewise, genome reduction of Mycobacterium lepraemurium from ancestral M. avium is likely unique to mice in causing murine leprosy (47).
Human M. leprae strains have been reported to infect a few animal species, such as nine-banded armadillos in the southern United States (8,48) and nonhuman primates (48,49). In the isolated infections of British red squirrels, the M. leprae strains were of human origin by hundreds of years while the M. lepromatosis strains bore a divergence of ;27,000 years from Mexican strain Mx1-22A (21). In red squirrels from continental northwest Europe, a recent study didn't find such infections (50). Therefore, considering the finding of M. lepromatosis in Asian patients with leprosy (17,18) and the likely introduction of the agent into Mexico 13,000 years ago during the earliest human Asia-Alaska-America migration (14), the most likely source of this bacillus in the British red squirrels could also be migrant humans from Asia ;27,000 years ago. In this regard, future studies on M. lepromatosis strains from Asian patients would be confirmative. Together, the new data and literature further corroborate our proposal.
The two leprosy bacilli share 71.4% of the genome but differ in 28.6%. The shared features should account for the dermal and neuronal invasionthe hallmark of leprosy. The preferential retention of PGL-I and its genes, favorable to bacterial survival, but inactivation of the detrimental PE/PPE components, are consistent features. Conversely, the variable genome features, in combination with the host genetic factors in diverse populations, should explain the variations in clinical and pathological manifestations of each infection, such as DLL and vascular occlusion associated with M. lepromatosis. Examples include three copies of the bifunctional diguanylate cyclase/phosphodiesterase gene in M. lepromatosis but one copy in M. leprae, 6 genes encoding the TetR transcriptional regulators in M. lepromatosis but 11 such genes in M. leprae, the presence of hemW in M. lepromatosis but absence in M. leprae, and so on. Clearly, the large numbers of shared and unique genome features require far more studies to delineate their relations with the infection. To this end, the annotated complete genome of M. lepromatosis has laid a foundation. In addition, recent report of passage of M. lepromatosis in mouse footpad (16) may provide viable bacteria for more studies as well as possible animal models in the near future.
The present study of the complete genome of M. lepromatosis also adds to the general knowledge on reductive genome evolution. This has been well documented in several bacterial endosymbionts and some other pathogens (51)(52)(53)(54), involving a shift in lifestyle from relative free-life to host restriction and intracellular living. These genomes undergo chromosomal instability, rearrangements, loss of dispensable genes through either pseudogenization by point mutations or large deletions, and a tendency for genome downsizing. In general, the rates of nucleotide substitution in CDS also increase drastically, except a few examples of low rates in some long-term endosymbionts (55)(56)(57)(58). Similar processes of massive pseudogenization have been detected in recent insect-associated bacterial symbionts due to their rich and stable new niches (59)(60)(61).

MATERIALS AND METHODS
Source of DNA and genome construction. The M. lepromatosis DNA from strain FJ924 was extracted from a smear slide that was prepared from autopsy liver tissue of a patient originally from Mexico (4). The The Complete Genome of M. lepromatosis Microbiology Spectrum smear had been dried, acid-fast stained (Kinyoun method) for microscopy, and archived for 6 years (2007-2013). Total genomic DNA was obtained with the QIAamp kit (Qiagen, Valencia, CA) following instructions from manufacturer. The extraction yielded ;3 ng DNA. A whole-genome library was then constructed using the KAPA kit (Kapa Biosystems, Wilmington, MA). The library was enriched by six PCR cycles, quantified, assessed for size distribution, and sequenced on the HiSeq 2000 sequencer with the 75-bp pair end configuration (Illumina, San Diego, CA). A total number of 69 million reads were obtained. The construction involved several steps. Human DNA contamination was removed first (14 million reads). From the remaining 55 million reads, mycobacterial DNA (11 million reads) were collected through BLAST matches (62) with M. leprae Br4923, filtering mostly contaminant DNA (44 million reads) from diverse other bacteria. The specific reads were assembled de novo using Velvet v1.2.10 (63), and the assembled contigs were aligned manually as well as using Bowtie v2.1.0 (64) to M. leprae template to build the first draft M. lepromatosis FJ924 genome. This draft was refined with GapFiller v1-10 (65) for gap closure and extension of contig edges, leading to announcement of a draft genome of 3,215,823 nucleotides with 39 contigs (66).
The published draft genome was used to map reads again, resulting in 12 million reads. These reads were assembled de novo with MIRA v2.1 with default parameters (67), and alignment of the assembled contigs led to third draft. The process of mapping, assembly, Gap filling, and draft construction was reiterated 5 times until seal of all blunt-end gaps. Eventually, a high-quality draft genome from 12.3 million reads and 525Â genome coverage, with 12 gaps that were all flanked by repetitive contigs, was obtained. Primers were designed at nonrepetitive regions of contig edges and PCRs along with Sanger sequencing of amplicons were used to close these gaps.
Genome annotation. The closed M. lepromatosis FJ924 genome was annotated independently with Basys (68), RAST web-serves (69), and Prokka v1.12 (70). Annotation results from the three methods were compared and mixed with BEACON v1.1 (71). Initial pseudogene prediction was based on Prokka annotation of adjacent CDS presenting the same annotation and a custom python script. In brief, the script uses LAST (72) alignments to detect CDS which are 75% shorter than a reference set of clustered proteomes at 95% amino acid identity (Table S8). Then, it mixes adjacent fragmented genes when they overlap with the same reference hit.
Pseudogenes annotations were revised exhaustively by using several types of BLAST analyses. The main one was based on BLASTX using the complete M. lepromatosis FJ924 genome sequence as a query against a protein database of M. haemophilum (E value 0.01). This revision produced length changes of several pseudogenes, fusion of others, and conversion of some initially annotated CDS to pseudogenes.
Comparative analysis of M. lepromatosis strains FJ924 and Mx1-22A. To examine the differences between the draft genome Mx1-22A and the complete genome of strain FJ924, the 126 contigs of Mx1-22A were queried to FJ924 in BLASTN searches (E value E-180). The results produced a complete alignment of 134 shared segments. The noncovered segments of FJ924 were identified and their annotation features were extracted (for those higher than 300 bp) (Table S1). SNPs between the two strains were identified with MAUVE (23). The initial SNP number of 314 was reduced to 160 after removing ambiguous nucleotides in the contigs of Mx1-22A. In the second round, all inaccuracies and SNPs within the end 500 nucleotides of the contigs were removed, leaving behind a set of 17 SNPs. These SNPs were checked, and six of them were removed due to their location in a short segment with many ambiguous nucleotides. A final set of 11 SNPs were studied.
Phylogenetics. Several available mycobacterial proteomes were selected for phylogenetic reconstruction. OrthoFinder v2. 3.3 (27) was used to detect, align, and concatenate 1,220 single core orthologous clusters of proteins shared between M. lepromatosis  (27,73,74). Obtained concatenated alignment was pruned with Gblocks v0.91b (75) with the option no gaps allowed. A total of 316,716 amino acid positions from 418,308 were maintained and used for phylogenetic reconstruction. The species tree was inferred by the Maximum Likelihood method (JTT1F1R3 model) using IQ-TREE v1.6.12 with 1000 ultrafast bootstrap and SH-aLRT as node support (74). An additional phylogenetic tree was inferred by the Maximum Likelihood method (JTT1G1I1F model and 500 bootstrap replicates) for 15 diguanylate cyclase/phosphodiesterase proteins of several mycobacterial species (601 sites) in MEGA7 (76).
A time tree was inferred by applying the RelTime method (77) using the previous concatenated alignment and phylogenetic tree (316,716 sites and 10 sequences), a model JTT1G1I1F and branch lengths estimation by the ordinary least-squares approach in MEGA11 (78). The calibration constraint was the minimum and maximum time of the node of divergence between the two leprosy bacilli (8.2-21.4 Mya) as previously reported (10).
Synteny and chromosomal rearrangement analysis. Alignment of several mycobacterial genomes was performed with MAUVE (23) using default parameters. GenoPlotR was used to plot Mauve results (79). The permutation file was converted with the program GRIMM in a distant matrix, which contained the minimal number of inversions required between a pair of genomes to explain differences in gene order (26). Distance matrix was charged in MEGA7 (76) and a neighbor joining algorithm was used to infer the rearrangement phylogeny (80). Orthology between CDS and pseudogenes from M. leprae Br4923 and M. lepromatosis FJ924 was obtained with a reciprocal BLASTN best hit strategy (E value = 1.0E-05). Nonreciprocal hits and those hits in disagreement with gene order were revised upon visual inspection. Synteny was the first criterion to assign orthology when more than one similar hit was detected allowing up to 5 hits in new BLASTN searches.
Identification and analysis of repetitive sequences. A BLASTN search using the genome of M. lepromatosis as both query and subject with an E value of 10 220 was performed to search for repetitive sequences. Segments with more than one overlapping hit were revised to produce the repeat segment with the largest length. Repetitive sequences of the same family were aligned in Ugene v33.0 (81). To determine the effect of repeat segments over chromosomal rearrangements, the 24 synteny breaks were delimited by the last CDS/pseudogene in a block and the first CDS/pseudogene in the next block with an orthologous in the genome of M. leprae. A plot comparing repetitive sequences and synteny breaks was made with Circos (82).
Identification of orthologous clusters. Orthologous clusters of proteins (genes) in the forms of core, pan, pairwise shared, and strain specific clusters were obtained with OrthoFinder v2.3.3 using the msa (mafft aligner) and IQ-TREE options (27,73,74). Data on the clusters were used for several analyses.
Estimation of the ages of pseudogenes. The relative ages of M. leprae and M. lepromatosis pseudogenes were estimated according to the method of Gomez-Valero et al. (7) with inclusion of the recent M. haemophilum genome data (12) for closer phylogeny. Briefly, this method used M. tuberculosis as outgroup to root the tree and compared M. haemophilum with either M. leprae or M. lepromatosis. The divergence between the two species determines the relative age of pseudogenization, ranging from 1, when it took place at the start of the divergence, to 0 at this time. The method is based on the estimation of the dN and dS (the number of nonsynonymous and synonymous substitutions per site, respectively) for each cluster of orthologous gene-pseudogene in the branches leading to M. haemophilum and a leprosy bacillus.
Estimation of the dN and dS values required use of clusters of orthologous genes. To produce clusters of orthology that incorporated M. leprae and/or M. lepromatosis pseudogenes and corresponding genes from both M. haemophilum and M. tuberculosis, a reciprocal best hit strategy was used with E value 1.0E-05. The set of pseudogenes/genes of M. lepromatosis and M. leprae were then compared, and they were further aligned to the sets of M. haemophilum and M. tuberculosis.
Alignments of genes and pseudogenes to maintain the codon structure, despite the indel mutations of the pseudogenes, were performed with MACSE (83). This program was also used to align clusters of functional genes in the four species. The program codeML (84) was used to obtain the dN and dS values as previously described (58). A custom python script was used to run those analyses.
To estimate relative pseudogene ages (Rpa), the following formulas were used, for M. lepromatosis The i (intercept) corresponds to the node of divergence between the leprosy bacilli and M. haemophilum lineages. The values dN ilp_ps and dN ile_ps are the numbers of nonsynonymous substitutions in the branches of a pseudogene from M. lepromatosis and M. leprae, respectively. To be pointed out, the term dN only shows a strict sense for the period of evolution as a gene.
The parameters f(dN ih ) and f'(dN ih ) are estimates, for orthologous genes, of the expected dN value for the branch of a leprosy bacillus as a function of the dN value for the branch of M. haemophilum considering that this gene had been evolving as a functional gene for the complete period of time. Orthologous genes evolve faster in the leprosy bacilli than in M. haemophilum, especially for dN (four times faster on average). By using dN values of 1,213 clusters of orthologous genes, the relationship between outcome and predictor, both logtransformed, was estimated by segmented linear regression. A two-piece linear relationship was assumed, namely, represented by two straight lines connected at an unknown breakpoint. Estimation was performed by means of Muggeo's algorithm (85)  The time of divergence in Mya between the lineage of leprosy bacilli and the M. haemophilum lineage (Rpa = 1) was estimated with the RelTime method (28). All mentioned scripts are available at https:// github.com/DiegoSantos-Garcia/Scripts. Data availability. This whole-genome shotgun project has been deposited at DDBJ/EMBL/GenBank under the accession number CP083405 and BioProject PRJNA281005.


Circadian Gene cry Controls Tumorigenesis through Modulation of Myc Accumulation in Glioblastoma Cells

Glioblastoma (GB) is the most frequent malignant brain tumor among adults and currently there is no effective treatment. This aggressive tumor grows fast and spreads through the brain causing death in 15 months. GB cells display a high mutation rate and generate a heterogeneous population of tumoral cells that are genetically distinct. Thus, the contribution of genes and signaling pathways relevant for GB progression is of great relevance. We used a Drosophila model of GB that reproduces the features of human GB and describe the upregulation of the circadian gene cry in GB patients and in a Drosophila GB model. We studied the contribution of cry to the expansion of GB cells and the neurodegeneration and premature death caused by GB, and we determined that cry is required for GB progression. Moreover, we determined that the PI3K pathway regulates cry expression in GB cells, and in turn, cry is necessary and sufficient to promote Myc accumulation in GB. These results contribute to understanding the mechanisms underlying GB malignancy and lethality, and describe a novel role of Cry in GB cells.

Introduction
Glioblastoma (GB) is the most common and aggressive type of glioma among all brain tumors, and it accounts for 57.3% of all gliomas [1]. It was classified in 2016 as a WHO grade IV diffuse oligodendroglial and astrocytic brain tumor, but the most recent classification (2021) includes these type of tumors in the "Gliomas, glioneuronal tumors, and neuronal tumors" group, termed as Glioblastoma, IDH-wildtype [2]. Despite current treatments, the median survival of GB patients is 15 months [3], and it is estimated that only 6.8% of patients survive five years after diagnosis [1]. To understand the genetic, molecular and cellular bases of gliomagenesis is fundamental for the development of effective therapies. In terms of histopathology and genetic expression, GB is a very heterogeneous type of tumor, even within the same patient [4]. However, there are common mutations in GB affecting different pathways that show mutual exclusion: the p53 pathway, the Rb pathway and components of the PI3K pathway [5].
Previous studies from our lab used a GB model in Drosophila, developed by Read and collaborators in 2009, that recapitulates key aspects of the disease both genetically and phenotypically [6][7][8][9][10][11][12][13]. This model is based on the expression of constitutively active forms of the epidermal growth factor receptor (EGFR λ ) and phosphatidyl inositol 3 kinase (PI3K) catalytic subunit (dp110 CAAX ) (orthologues of EGFR and PI3K catalytic subunit in Drosophila, respectively). We used the binary expression system Gal4/UAS [14] to express EGFR λ and PI3K dp110 CAAX specifically in glial cells under the control of repo-Gal4 driver [6]. The co-activation of EGFR and PI3K signaling pathways in Drosophila glial cells reproduces the cascade of signaling events that occurs in GB patients [6]. In consequence, GB cells upregulate myc expression, which is essential for tumoral transformation, and the glial tumor cell numbers increase along with the expansion of the glial membrane. As a result, Int. J. Mol. Sci. 2022, 23, 2043 2 of 16 GB progression causes a reduction in the number of synapses in neighboring neurons and premature death [6,9,15]. Furthermore, EGFR and PI3K pathway co-activation regulates processes such as progression and entry into the cell cycle and protein synthesis [6,7].
c-myc is one of the oncogenes most amplified in human cancer, including GB. About 60%-80% of human GB cases show elevated Myc levels [16]. Myc regulates cell proliferation, transcription, differentiation, apoptosis and cell migration. It is the point where EGFR and PI3K pathways converge; thus, Myc is considered essential for GB transformation [6,[16][17][18]. Furthermore, in vitro and in vivo studies have shown that myc inhibition prevents glioma formation, inhibits cell proliferation and survival and even induces disease regression [16,19]. These features are conserved in Drosophila [6].
In the recent years, the study of alterations in circadian rhythm genes has emerged in different types of cancer, including GB [20]. Previous reports suggested that circadian rhythm genes play essential roles in different aspects of tumor progression. The central clock organizes the oscillations and rhythmicity of the physiological processes and modulates the expression of genes related to cell proliferation or differentiation, such as cell cycle components [21], proto-oncogenes and tumor suppressors [22].
In mammals, the structure responsible for coordinating circadian behavior throughout the body is the suprachiasmatic nucleus (SCN), located in the anterior region of the hypothalamus and made up of about 50,000 neurons in humans [23]. All the neurons that compose the central clock express the core circadian genes that control the oscillations that organize the cycles of the whole organism in absence of environmental cues. Furthermore, synchronization of the internal clock with light/dark cycles relies on cryptochrome protein (Cry), a blue light photopigment expressed in certain subsets of clock neurons. Cry is a receptor of near-UV/blue light and a regulator of gene expression that belongs to the group of DNA photolyases. It was suggested that the last universal common ancestor (LUCA) had one or several photolyases, supporting the evolutionary conservation of cryptochrome genes [24]. However, the mammalian gene that plays the role of Drosophila cry remains unknown. Interestingly, Drosophila Cry also acts as the mammalian Cry when expressed in peripheral clocks [25]. Besides, cry1 expression is androgen responsive, Cry1 regulates DNA repair and the G2/M transition and it is associated with poor outcome in prostate cancer and colorectal cancer.
Regarding GB, studies in patients with primary gliomas show an association between a specific per1 variant with overall glioma risk. Several circadian genes, including cry1, exhibited differential expression in GB samples compared to control brains as described in the literature [26,27], and in human cancer gene expression databases (https://www.proteinatlas. org, accessed on 1 February 2022; https://cancer.sanger.ac.uk/, accessed on 10 January 2022). Besides, the expression of the circadian gene clk is significantly enhanced in high-grade gliomas and correlates with tumor progression [28]. Moreover, per1 and per2 expression increases the efficacy of radiotherapy also in GB cells [29].
Furthermore, high levels of cry1 inversely correlate with median survival in GB patients, acting as signal of poor prognosis (http://gepia.cancer-pku.cn/detail.php?gene= CRY1, accessed on 1 February 2022). Still, the functional mechanism of Cry in cancer susceptibility and carcinogenesis remains unsolved.
Different studies show a relationship between Cry and Myc [30]; c-Myc levels decrease in cry1/cry2 null mutant mice [31]. Besides, cry1 expression is induced by Myc in GB cells in culture [32].
Taking into account the deregulation in the expression of circadian genes in tumor tissues and the pre-established relationship between Cry and myc, which is a key player in GB, here we show that cry is regulated by PI3K pathway, cry expression enhances Myc accumulation in GB cells and it is necessary for GB progression.

Cry Expression in Glioblastoma
To determine if cry expression was affected in glioma samples, we extracted RNA from the heads of 7-day-old adult control and glioma flies. Quantitative RT-PCR results (see Table 1 in Materials and Methods) indicate that cry mRNA levels are 50 times higher in glioma samples as compared to controls ( Figure 1A). This result goes in line with the data retrieved from TCGA-GBM dataset (at http://gliovis.bioinfo.cnio.es/, accessed on 1 February 2022) that indicate a significant increase of cry1 mRNA levels (RNA-seq) in GB samples, as compared to non-tumor tissue.
Next, to determine if cry upregulation occurs in GB cells, we used a specific reporter line that generates a green fluorescent protein (GFP) tagged form of Cry (GFP-Cry) and visualized adult brains in confocal microscopy. The images show the GFP signal (Cry) and glial membrane marked in red with myristoylated red fluorescent protein (mRFP) ( Figure 1B-E,B'-E'). The quantification of GFP-Cry and mRFP co-localization is higher in glioma samples than in controls ( Figure 1B,C,F) suggesting an accumulation of Cry in glioma cells. This signal is restored to control levels upon cry knockdown by means of RNAi expression in glial or glioma cells ( Figure 1D-F).
Next, we analyzed human mRNA expression databases for Glioblastoma multiforme (http://gliovis.bioinfo.cnio.es/, accessed on 10 January 2022). The results indicate that cry in GB patients is transcriptionally upregulated in primary tumors ( Figure 1G) and cry1 upregulation correlates with worse prognosis ( Figure 1H). Moreover, cry1 is also upregulated in secondary GB ( Figure 1I) and correlates with poor prognosis in secondary GB patients ( Figure 1J). All together, these results indicate that cry is transcriptionally upregulated in GB cells in Drosophila and patients and suggest a role in GB malignancy and aggressiveness.

Cry Mediates GB Progression and Neurodegeneration
To determine the contribution of cry to GB progression, we used a previously validated protocol to quantify tumor progression and the associated neurodegeneration in Drosophila [7, 9,11]. We stained adult control brains and compared them with GB, GB + cryRNAi and wt brains expressing cryRNAi in glial cells. We used a specific antibody against repo to visualize the nuclei of all glial cells and quantified the number of glial cells in the confocal images (Figure 2A-E). The results indicate that GB samples have a significant increase in the number of glial cells compared to control samples, but this increase depends on cry expression (Figure 2A-C,E). Besides, knockdown of cry in normal glia does not alter the number of glial cells ( Figure 2D,E). In addition, we quantified the volume of glial membrane. We used Imaris software to measure the volume of the red signal that corresponds to a myristoilated form of RFP (mRFP) expressed in glial cells under the control of repo-Gal4. The quantification of the volume show a significant expansion of glial membrane in GB compared to control samples, but this increase depends on cry expression ( Figure 2A'-C',F). Again, knockdown of cry in normal glia does not alter the volume of glial membrane ( Figure 2D',F). These results suggest that cry expression is required for GB progression, but not for normal glia development.

cry Expression in Glioblastoma
To determine if cry expression was affected in glioma samples, we extracted RNA from the heads of 7-day-old adult control and glioma flies. Quantitative RT-PCR results (see Table 1 in Materials and Methods) indicate that cry mRNA levels are 50 times higher in glioma samples as compared to controls ( Figure 1A). This result goes in line with the data retrieved from TCGA-GBM dataset (at http://gliovis.bioinfo.cnio.es/, accessed on February 1 st 2022) that indicate a significant increase of cry1 mRNA levels (RNA-seq) in GB samples, as compared to non-tumor tissue. Next, we studied the impact of GB progression and cry expression in neighboring neurons. We counted the number of synapses in motor neurons of adult neuromuscular junction (NMJ), a standardized tissue to study neurodegeneration [9,11,33]. To visualize synapses, we used an anti-Brp antibody (nc82) to detect active zones in the neurons, and counted the number of synapses in control samples, GB, GB + cryRNAi and normal glia + cryRNAi ( Figure 2G-J). The quantification of synapse number ( Figure 2K) shows that GB induction provokes a significant reduction in the number of synapses as compared to control samples, compatible with a neurodegenerative process. This effect was previously described [7, 9,11] as a consequence of GB progression. Moreover, cry knockdown in GB prevents the reduction in the number of synapses, and cryRNAi expression in normal glial cells does not cause any detectable change in the number of synapses. Finally, we aimed to determine the systemic effect of cry. We expressed cryRNAi in glia or GB cells, and we analyzed the life span of adult flies. The results show that GB causes a significant reduction of life span and a premature death, which is prevented by cryRNAi expression in GB cells. Moreover, cryRNAi expression in normal glial cells does not reduce lifespan but causes a significant increase in the average lifespan ( Figure 2L).

Signaling Pathway to Control Cry Upregulation
To decipher the specific signaling pathway responsible for cry transcriptional activation in GB cells, we analyzed the contribution of the two main pathways activated in this model of GB, EGFR and PI3K. Both pathways converge in the expression of the gene myc (see Figure 3A for detailed genetic epistasis in GB). Thus, we analyzed the contribution of PI3K, EGFR and myc to cry upregulation. We measured the fluorescent signal of GFP-cry reporter in control adult brains ( Figure 3B-B") and compared it with adult brains upon expression of the constitutively active forms of PI3K ( Figure 3C-C") or EGFR ( Figure 3D-D") in glial cells (under the control of repo-Gal4). In addition, we analyzed the GFP-cry signal in glial cells upon myc upregulation ( Figure 3E-E"). We quantified in the confocal images the signal of GFP that co-localizes with glial membranes (mRFP) ( Figure 3F). The results indicate that PI3K expression is sufficient to increase GFP-cry signal but not EGFR or myc overexpression. These results suggest that PI3K upregulation induces cry transcription, and EGFR or myc expression do not induce cry expression in glial cells.

Cry Regulates Myc Expression in Glial Cells
Next, to determine the epistatic relation between cry and myc, we analyzed Myc protein accumulation in glial cells upon cry expression. First, to analyze if Cry is sufficient to cause an increase in Myc protein levels, we used a specific antibody against Myc and analyzed Myc signal levels upon cry overexpression, myc overexpression or cry + myc overexpression in glial cells ( Figure 4A-D'). The quantification of Myc surface signal that coincides with glial cells (anti-repo) showed that cry expression in glia is sufficient to increase Myc protein signal in glial cells, comparable to myc upregulation. In addition, cry + myc upregulation show a summation effect on the increase of Myc protein levels ( Figure 4E). To conclude if cry is required for myc expression in GB, we quantified glial Myc signal in the control, cry RNAi, GB, GB + cryRNAi and cry upregulation ( Figure 4F-J'). The quantifications indicate that cryRNAi in glial cells does not reduce the amount of Myc in glial cells ( Figure 4K). In addition, GB condition triggers the number of Myc positive glial cells, as well as cry upregulation in glial cells ( Figure 4K). Finally, cryRNAi expression in GB cells prevents the accumulation of Myc in GB cells. Taking all these results together, we conclude that cry is sufficient to trigger Myc accumulation in glial cells, and cry expression is necessary for Myc accumulation in GB condition.

cry Regulates myc Expression in Glial Cells
Next, to determine the epistatic relation between cry and myc, we analyzed Myc protein accumulation in glial cells upon cry expression. First, to analyze if Cry is sufficient to cause an increase in Myc protein levels, we used a specific antibody against Myc and analyzed Myc signal levels upon cry overexpression, myc overexpression or cry + myc overexpression in glial cells ( Figure 4A-D'). The quantification of Myc surface signal that coincides with glial cells (anti-repo) showed that cry expression in glia is sufficient to increase Myc protein signal in glial cells, comparable to myc upregulation. In addition, cry + myc upregulation show a summation effect on the increase of Myc protein levels ( Figure  4E). To conclude if cry is required for myc expression in GB, we quantified glial Myc signal in the control, cry RNAi, GB, GB + cryRNAi and cry upregulation ( Figure 4F-J'). The quantifications indicate that cryRNAi in glial cells does not reduce the amount of Myc in glial cells ( Figure 4K). In addition, GB condition triggers the number of Myc positive glial cells,

Cry Contribution to GB Progression
To investigate the contribution of Cry to glioma progression, we determined the number of glial cells and volume of glial membrane network in control adult brain, GB (PI3K + EGFR), PI3K + cry, EGFR + cry or myc + cry expressed in glial cells ( Figure 5A-E'). The quantification showed that all these genetic combinations cause an increase in the number of glial cells as compared to control brains ( Figure 5F). However, only the GB condition provoked an expansion of the glial membrane volume, and the combination of PI3K + cry, EGFR + cry or myc + cry showed a volume of glial membrane comparable to control brains ( Figure 5G). To further determine the contribution of cry to GB expansion, we analyzed the contribution of single gene upregulation in glial cells for cry or myc, and the combination of cry + myc expression ( Figure 5H-K'). The quantification of glial cell number showed that cry or myc expression alone, or in combination, is sufficient to increase the number of glial cells with respect to control samples ( Figure 5L). Nevertheless, none of these genetic modifications is sufficient to expand glial membrane volume ( Figure 5M).
These results suggest that cry or myc are sufficient to trigger glial cell number increase in adult brains, but not to expand the volume of glial membrane network. as well as cry upregulation in glial cells ( Figure 4K). Finally, cryRNAi expression in GB cells prevents the accumulation of Myc in GB cells. Taking all these results together, we conclude that cry is sufficient to trigger Myc accumulation in glial cells, and cry expression is necessary for Myc accumulation in GB condition.

cry Contribution to GB Progression
To investigate the contribution of Cry to glioma progression, we determined the number of glial cells and volume of glial membrane network in control adult brain, GB (PI3K + EGFR), PI3K + cry, EGFR + cry or myc + cry expressed in glial cells ( Figure 5A-E'). The quantification showed that all these genetic combinations cause an increase in the number of glial cells as compared to control brains ( Figure 5F). However, only the GB condition provoked an expansion of the glial membrane volume, and the combination of PI3K + cry, EGFR + cry or myc + cry showed a volume of glial membrane comparable to control brains ( Figure 5G). To further determine the contribution of cry to GB expansion, we analyzed the contribution of single gene upregulation in glial cells for cry or myc, and the combination of cry + myc expression ( Figure 5H-K'). The quantification of glial cell number showed that cry or myc expression alone, or in combination, is sufficient to increase the number of glial cells with respect to control samples ( Figure 5L). Nevertheless, none of these genetic modifications is sufficient to expand glial membrane volume ( Figure  5M). These results suggest that cry or myc are sufficient to trigger glial cell number increase in adult brains, but not to expand the volume of glial membrane network.

Cry Upregulation in Glial Cells Causes Synapse Loss and Premature Death
It was previously described that GB progression induces synapse loss, an early symptom of neurodegeneration. To determine the contribution of cry to synapse loss, we counted the number of active zones in motor neurons of adult neuromuscular junction in the control, GB (PI3K + EGFR), PI3K + cry, EGFR + cry or myc + cry samples ( Figure 5N-R). The quantification of the number of active zones showed that the expression in glial cells of GB (PI3K + EGFR), PI3K + cry, EGFR + cry or myc + cry is sufficient to reduce the number of synapses in NMJ neurons ( Figure 5S).
Finally, to evaluate the systemic effect of GB and glial expression of PI3K + cry, EGFR + cry or my c + cry, we analyzed the lifespan of adult individuals. The results show that GB causes a premature death, as previously described in Drosophila and mice Xenografts [8,9,11], glial upregulation of EGFR + cry or myc + cry causes a significant reduction of lifespan but less aggressive than GB, and PI3K + cry upregulation in glial cells does not reduce lifespan ( Figure 5T).

Discussion
Different studies have established a relation between alterations in circadian rhythm genes and cancer [32,34]. Specifically, one of the genes associated with different types of cancer is cry [35][36][37]. Thus, this study aims to investigate the role of cry in a Drosophila GB model.
The previous work of Luo et al. 2012 [38] describes a reduction of the number of glial cells positive for cry1/2 expression in glioma tissue compared to normal tissue. However, the authors show that glioma cells that are positive for Cry1/2 show an increase in the amount of Cry1/2 with respect to non-tumoral tissue. Moreover, both Madden et al. 2014 [26] (with a sample 10 times larger than that of Luo et al. 2012) and Wang et al. 2021 [27] (using data from three different databases) analyzed the expression of circadian genes in glioma tissue compared to healthy tissue and conclude that cry1 is overexpressed in glioma tissue. We also found this result in other databases such as https://www.proteinatlas.org/ and http://gliovis.bioinfo.cnio.es/, accessed on 10 January 2022, [39], which in turn is compatible with the observations in the Drosophila model of GB.
Nonetheless, Fan et al. [40] investigated the role of Cry2 in rat glioma cells and observed that cry2 mRNA and protein levels showed aberrant rhythmic periodicity of 8 h, compared to 24 h in normal tissue. Thus, futures studies on the contribution of circadian rhythms genes should take into consideration the variations of expression.
On the contrary, Dong et al. [41] state that glioblastoma stem cells (GSCs) displayed robust circadian rhythms dependent on core clock transcription factors. The use of Cry1/2 agonists induced anti-tumor effects suggesting that GSCs are sensitive to cry1/2 activity. Taken the different conclusions into consideration, most of the literature and our data suggest that cry is upregulated in glioma cells and promotes glioma progression; however, the role of cry expression in Drosophila, or cry1/2 expression in mammals, may differ according to the glioma subtype, the specific mutations in glioma cells and the cell population of study within the glioma and the hour of the day.
We described an increase in cry1 mRNA levels in human GB samples and in a well-studied Drosophila model of GB. However, we cannot conclude that GB cells show an increase in cry transcription, or an enhancement of cry mRNA stability. The Drosophila GB model is based on the activation of the two most frequently mutated pathways in GB, PI3K and EGFR, which converge in Myc as a coincidence point. These pathways are of great relevance to promote GB cells expansion, GB progression and, in consequence, the deterioration of neighboring neurons and a premature death. The results indicate that cry upregulation in Drosophila GB cells depends on PI3K expression, and it is required for GB cells number increase and synapse loss ( Figure 6). In addition, cry expression in glial cells is sufficient to increase the number of glial cells. However, cry expression is expendable for normal wt glial growth during development. Taking into consideration that cry is upregulated in GB cells and promotes glia cells number increase, we did not observe any contribution to normal glia development, which makes Cry a potential target for GB treatment.
Besides, we show that Cry is necessary and sufficient to induce myc expression in GB cells. This agrees with in vitro studies that revealed an increase in Myc levels as a result of cry upregulation [32]. Therefore, we propose that cry is part of the PI3K-Myc signaling pathway in GB, where cry upregulation would be associated with glial cells number increase. However, PI3K is a highly promiscuous enzyme that participates in numerous signaling pathways, and the results suggest that Cry contribution is restricted to the malignant features of GB dependent on myc, such as GB cell number increase and neurodegeneration. However, cry expression is independent of glial membrane expansion characteristic of GB progression. Besides, cry expression in glial cells partially reduces lifespan, but is less aggressive than GB. This result suggests that Cry plays a central role in GB and is required for GB formation, and cry mutations might be responsible for several features of GB. The human gene expression databases indicate that cry1 expression levels correlate negatively with lifespan, and it is associated with a poor prognosis. In conclusion, these results suggest that further studies on the contribution of Cry1 to human GB progression could lead to novel strategies to treat GB patients. great relevance to promote GB cells expansion, GB progression and, in consequence, the deterioration of neighboring neurons and a premature death. The results indicate that cry upregulation in Drosophila GB cells depends on PI3K expression, and it is required for GB cells number increase and synapse loss ( Figure 6). In addition, cry expression in glial cells is sufficient to increase the number of glial cells. However, cry expression is expendable for normal wt glial growth during development. Taking into consideration that cry is upregulated in GB cells and promotes glia cells number increase, we did not observe any contribution to normal glia development, which makes Cry a potential target for GB treatment. Besides, we show that Cry is necessary and sufficient to induce myc expression in GB cells. This agrees with in vitro studies that revealed an increase in Myc levels as a result Recent publications describe the communication between GB cells and neurons in human GB cells and mice xenografts based on the establishment of electrical and chemical synapses, which are essential for tumor progression [42,43]. A possible explanation for GB prevention by cry downregulation arises from the non-circadian function of Cry as a regulator of synapse number through the genetic and physical association with the key presynaptic protein Bruchpilot (Brp) [44,45]. In Drosophila, cry mutants show reduced brp expression levels. Actually, Cry interacts physically with Brp to modulate its stability, and triggers its degradation activated by light. Therefore, it is possible that cry overexpression in glial cells promotes the establishment of synapses with neurons. Moreover, the absence of light input impairs Brp degradation in glial cells, thus promoting tumoral progression.
In conclusion, further experiments are required to unveil the molecular interactions of Cry and Brp proteins, including the putative formation of abnormal synapses between glial and neurons under the GB condition.
The studies of other groups describe the beneficial effects of haloperidol on cry1 expression in GB cells, but these results obtained in cell culture suggest that the doses required to treat patients might be toxic; in consequence, specific delivery strategies combined with haloperidol are worth of study. In addition, we observed significant effects of cry knockdown in normal glial cells, in line with Bolukbasiet al., who recently described the extension of lifespan by foxo upregulation in glial cells [46]. We observed an effect of cry upregulation in the number of glial cells ( Figure 5L). Given that cry and foxo respond to PI3K pathway, it is tempting to speculate that cry expression is relevant for lifespan extension by PI3K pathway, and associated behaviors such as diet restriction.
The classical definition of Cry as a regulator of circadian rhythms can now be expanded to the biology of glial cells, GB progression and the expansion of lifespan in Drosophila. This plethora of different phenotypes associated with one gene is now a common feature previously described for Troponin I [47][48][49], Caspases [50][51][52] or even other circadian genes as per1 [53,54] and contributes to the explanation of the multiple phenotypes observed in patients.
This study describes the epistatic relationship between PI3K, cry and myc and the relevance for GB progression. The strengths of this study rely on the importance of understanding the mechanisms underlying the progression of a fatal tumor as GB, and the reliability of Drosophila as an animal model useful to study human disease. However, it is important to take into consideration the limitations of the results to put them in perspective. We used a model based on the activation of PI3K and EGFR pathways that reproduces the key features of human disease progression, but the contribution of additional mutations such as IDH or TP53 in GB require further studies; thus, new models in flies or other animal models of study will contribute to validate and narrow down our findings.
The stock containing UAS-cryRNAi was previously generated and validated [57]; this construct produces a double-stranded RNA that corresponds to the 300-799 region of cryRA mRNA.
The glioma-inducing line contains the UAS-dEGFR λ and UAS-dp110 CAAX transgenes that encode for the constitutively active forms of the human orthologues PI3K and EGFR, respectively [6]. The Repo-Gal4 line drives the Gal4 expression to glial cells and precursors [59,60] combined with the UAS-dEGFR λ and UAS-dp110 CAAX line allow us to generate a glioma thanks to the Gal4 system [14]. To visualize glial or GB cells membrane, we induced the expression of a myristoylated form of red fluorescent protein (UAS-mRFP, described in [9]) under the control of the specific glial promoter repo-Gal4.
Gal80 TS is a repressor of the Gal4 activity at 18 • C, although at 29 • C is inactivated [61]. The tub-Gal80 TS construct was used in all the crosses to avoid the lethality caused by the glioma development during the larval stage. The crosses were kept at 17 • C until the adult flies emerged. To inactivate the Gal80 TS protein and activate the Gal4/UAS system to allow for the expression of our genes of interest; the adult flies were maintained at 29 • C for 7 days except in the survival assay (flies were at 29 • C until death).

Immunostaining and Image Acquisition
All tissues were treated in simultaneously for each experiment. Adult brains were dissected and fixed with 4% formaldehyde in phosphate-buffered saline for 20 min, whereas adult NMJ were fixed for 10 min; in both cases, samples were washed 3 × 15 min with PBS + 0.4% triton, blocked for 1 h with PBS + 0.4% triton + BSA 5%, incubated overnight with primary antibodies, washed 3 × 15 min, incubated with secondary antibodies for 2 h and mounted in Vectashield mounting medium with DAPI in the case of the brains. The primary antibodies used were anti-Repo mouse (1/200; DSHB, Iowa City, IA, USA) to recognize glial nuclei, anti-Bruchpilot-nc82-mouse (1/50; DSHB, Iowa City, IA, USA) to recognize the presynaptic protein Bruchpilot, anti-HRP rabbit (1/400; Cell Signaling, Danvers, MA, USA) to recognize neuronal membranes, anti-GFP rabbit (1:500; DSHB, Iowa City, IA, USA) and anti-Myc guinea pig (1/100; DSHB, Iowa City, IA, USA) to recognize the nuclear protein Myc. The secondary antibodies used were anti-mouse, -rabbit or -guinea pig Alexa 488 or 647 (1/500; Life Technologies, Carlsbad, CA, USA). Images were taken by a Leica SP5 confocal microscopy applying same conditions for each experiment.

qRT-PCR
The mRNA for all samples was extracted from adult brains and processed in parallel. For this, 1-to 4-day-old male adult mice were maintained at 29 • C for 7 days and collected on dry ice at ZT6. Total RNA was extracted by triplicate from 30 heads. RNA was extracted with TRIzol and phenol chloroform. cDNA was synthetized from 1 µg of RNA and cDNA samples from 1:5 dilutions were used for real-time PCR reactions. Transcription levels were determined in a 14 µL volume in duplicate using SYBR Green (Applied Biosystems, Waltham, MA, USA) and 7500 qPCR (Thermo Fisher Scientific, Waltham, MA, USA). We analyzed transcription levels of cry using Rp49 as a housekeeping gene reference.
Sequences of primers were as follows. After completing each real-time PCR run, with cycling conditions of 95 • C for 10 min, 40 cycles of 95 • C for 15 s and 55 • C for 1 min, and outlier data were analyzed using 7500 software (Applied Biosystems, Waltham, MA, USA). Ct values by triplicate of duplicates from three biological samples were analyzed calculating 2DDCt.

Survival Assays
Lifespan was determined under 12:12 h LD cycles at 29 • C conditions. Three replicates of 30 1-to 4-day-old male adults were collected in vials containing standard Drosophila media and transferred every 2-3 days to fresh Drosophila media.

Quantification
Fluorescent reporter-relative cry signals within brains were determined from images taken at the same confocal settings avoiding saturation. For the analysis of co-localization rates, "co-localization" tool from LAS AF Lite software (Leica, Wetzlar, Germany), was used taking the co-localization rate data for the statistics analyzing the co-localization between the green signal (both cases) and signal coming from glial tissue from three slices per brain in similar positions of the z axis.
Glial network was marked by a UAS-myristoylated-RFP reporter (mRFP) specifically expressed under the control of repo-Gal4. The total volume was quantified using the Imaris surface tool (Imaris 6.3.1 software, Oxford Instruments, Abingdon, UK). Glial nuclei were marked by staining with the anti-Repo (DSHB). The number of Repo + cells and number of synapses (anti-nc82; DSHB) were quantified by using the spots tool in Imaris 6.3.1 software (Oxford Instruments, Abingdon, UK). We selected a minimum size and threshold for the spot in the control samples of each experiment: 0.5 µm for active zones and 2 µm for glial cell nuclei. Myc glial signal was quantified using the Imaris surface tool (Imaris 6.3.1 software, Oxford Instruments, Abingdon, UK) creating a mask for the glial nuclei signal and exclusively selecting the myc signal corresponding to glial nuclei. Then we applied the same conditions to the analysis of the corresponding experimental sample.

Statistics
The results were analyzed using the GraphPad Prism 5 software. Quantitative parameters were divided into parametric and nonparametric using the D'Agostino and Pearson omnibus normality test, and the variances were analyzed with F test. The t-test and ANOVA test with Bonferroni's post hoc were used in parametric parameters, using Welch's correction when necessary. The survival assays were analyzed with Mantel-Cox test. The p limit value for rejecting the null hypothesis and considering the differences between cases as statistically significant was p < 0.05 (*). Other p-values are indicated as ** when p < 0.01 and *** when p < 0.001.

Human GB Databases
We used a public open access database (http://gliovis.bioinfo.cnio.es/, accessed on 1 February 2022) to analyze the expression of human Cry1 gene in GB samples. We used the "Adult" samples in CGGA Dataset and included the data from primary and secondary tumor types. The data shown in Figure 1 correspond to the "expression" and "survival" tabs. Please note that nomenclature corresponds to the 2016 classification. GBM-Glioblastoma multiforme.