Maryanna Cristiano Simão Exaptação de elementos de transposição (TEs) em Drosophila mojavensis e viés de expressão de genes e TEs em gônadas de D. mojavensis e D. arizonae e seus híbridos São José do Rio Preto 2022 Câmpus de São José do Rio Preto Maryanna Cristiano Simão Exaptação de elementos de transposição (TEs) em Drosophila mojavensis e viés de expressão de genes e TEs em gônadas de D. mojavensis e D. arizonae e seus híbridos São José do Rio Preto 2022 Tese apresentada como parte dos requisitos para a obtenção do título de Doutor em Biociências, junto ao Programa de Pós-Graduação em Biociências, Área de Genética e Biologia Evolutiva, do Instituto de Biociências, Letras e Ciências Exatas da Universidade Estadual Paulista “Júlio de Mesquita Filho”, Câmpus de São José do Rio Preto. Financiadores: CAPES IDEX LYON Orientadora: Profª. Drª. Claudia Marcia Aparecida Carareto Coorientadora: Profª. Drª. Cristina Vieira N°d’ordre NNT : 2021LYSE1327 THESE de DOCTORAT DE L’UNIVERSITE DE LYON Opérée au sein de L’Université Claude Bernard Lyon 1 Ecole Doctorale N° ED 341 - Évolution, Écosystèmes, Microbiologie, Modélisation (E2M2) Spécialité de doctorat : Discipline : (Eventuellement) Soutenue publiquement/à São José do Rio Preto, São Paulo, Brésil le 17/12/2021, par : Maryanna Cristiano Simão Exaptation d’éléments transposables (TEs) chez Drosophila mojavensis et biais d'expression des gènes et des TEs chez les gonades de D. mojavensis et D. arizonae et leurs hybrides _________________________________________________________________ Devant le jury composé de : Lilian Ravazzi, Prof, UNESP (qui pourrait être présidente), Université d’état Paulista São José do Rio Preto (Brésil) : Examinatrice Maria del Pilar Garcia Guerreiro, Professeure Université Autonome de Barcelone (Espagne) : Rapporteur Elgion Lucio da Silva Loreto, Professeur, Université Fédérale de Santa Maria (Brésil) : Rapporteur/Examinateur Gustavo Campos e Silva Kuhn, Professeur, Université du Minas Gerais (Brésil) : Rapporteur/Examinateur Marie Fablet, MCU HDR, Université Lyon1 (France): Examinatrice Claudia Marcia Aparecida Carareto Professeur Université d’état Paulista São José do Rio Preto : Co-Directrice de thèse Cristina Vieira, Professeure, Université Lyon1 (France) : Directrice de thèse Ficha catalografica gerada pela autora Essa ficha não pode ser modificada Simão, Maryanna Cristiano Exaptação de elementos de transposição (TEs) em Drosophila mojavensis e viés de expressão de genes e TEs em gônadas de D. mojavensis e D. arizonae e seus híbridos / Exaptation d’éléments transposables (TEs) chez Drosophila mojavensis et biais d'expression des gènes et des TEs chez les gonades de D. mojavensis et D. arizonae et leurs hybrides / Maryanna Cristiano Simão. – São José do Rio Preto, 2022 139 p. : il., tabs., maps Tese (doutorado) – Universidade Estadual Paulista (UNESP), Instituto de Biociências, Letras e Ciências Exatas, São José do Rio Preto Orientadora: Claudia Marcia Aparecida Carareto Tese (doutorado) – Université Claude Bernard Lyon 1, Lyon Orientadora: Cristina Vieira 1. Elementos de Transposição. 2. Exaptação. 3. Genes com viés sexual. 4. Híbridos. 5. Drosophila. I. Título. Maryanna Cristiano Simão Exaptação de elementos de transposição (TEs) em Drosophila mojavensis e viés de expressão de genes e TEs em gônadas de D. mojavensis e D. arizonae e seus híbridos Comissão Examinadora Profª Drª Claudia Marcia Aparecida CARARETO UNESP – Câmpus de São José do Rio Preto – Orientadora Profª Drª Cristina VIEIRA Université Claude Bernard Lyon1 – Orientadora Profª Drª Lilian MADI-RAVAZZI UNESP – Câmpus de São José do Rio Preto Prof. Dr. Elgion Lucio da Silva LORETO Universidade Federal de Santa Maria (UFSM) Prof. D.r Gustavo Campos e Silva KUHN Universidade Federal de Minas Gerais (UFMG) MCU HDR Drª Marie FABLET Université Claude Bernard Lyon1 São José do Rio Preto 17 de dezembro de 2021 Tese apresentada como parte dos requisitos para a obtenção do título de Doutor em Biociências, junto ao Programa de Pós-Graduação em Biociências, Área de Genética e Biologia Evolutiva, do Instituto de Biociências, Letras e Ciências Exatas da Universidade Estadual Paulista “Júlio de Mesquita Filho”, Câmpus de São José do Rio Preto. Financiadores: CAPES IDEX LYON Orientadora: Profª. Drª. Claudia Marcia Aparecida Carareto Coorientadora: Profª. Drª. Cristina Vieira Este trabalho foi realizado sob convenção de co-tutela entre Universidade Estadual Paulista (UNESP) - Brasil e l’Université Claude Bernard (Lyon1 - UCBL) – França, no Laboratório de Evolução Molecular, do Instituto de Biociências, Letras e Ciências Exatas (IBILCE/UNESP) –São José do Rio Preto/SP, e no Laboratoire de Biométrie et Biologie Évolutive (LBBE/Lyon1). Universidade Estadual Paulista (UNESP – IBILCE) Laboratório de Evolução Molecular Departamento de Biologia Rua Cristóvão Colombo, 2265, Jardim Nazareth 15054-000 São José do Rio Preto-SP Brasil Université Claude Bernard – Lyon 1 Laboratoire de Biométrie et Biologie Evolutive CNRS UMR 5558 43 Boulevard du 11 novembre 1918 69622 Cedex Villeurbanne Maryanna CRISTIANO SIMÃO maryanna.simao@unesp.br simaomc.bio@gmail.com maryanna.simao@etu.univ-lyon1.fr Palavras chaves – Mots clés – Keywords Elementos de Transposição - Eléments transposables – Transposable elements Exaptação – Exaptation - Exaptation Domínios zf-BED - Domaines zf-BED, zf-BED domains Gene lncRNAs - Gènes lncRNAs – lncRNAs genes Genes com viés sexual - Gènes liés au sexe -Sex-biased genes Híbridos - Hybrides – Hybrids Drosophila mojavensis Drosophila arizonae À Maria, minha avó, melhor amiga, e maior saudade, que me ensinou o significado de força. AGRADECIMENTOS Aos meus pais, Sylas Simão dos Santos e Lucelaine Aparecida Cristiano Simão dos Santos, e minha irmã Maraynna Cristiano Simão por estarem ao meu lado, tentando entender e me apoiar durante alguns momentos mais difíceis. Às minhas orientadoras, Claudia Marcia Aparecida Carareto e Cristina Vieira, pela paciência, insistência e perseverança comigo durante todo esse processo. Obrigada por todo o tempo compartilhado, pela Ciência, pelas conversas e por mostrar a diferença entre o caminho mais fácil e o caminho que se deve ser feito. Obrigada por permitirem eu fazer parte desse grupo de mulheres que fazem Ciência e que mudam pessoas. Aos amigos que eu fiz ao longo dos anos no IBILCE/UNESP, sem os quais não seria possível a finalização dessa jornada, dentre os quais destaco minha amiga e irmã Geisla de Oliveira por conhecer todos meus defeitos continuar ao meu lado. Meu companheiro e parceiro Gabriel Spanghero que participou dos meus melhores e piores devaneios e reverberações ao longo desses anos. Aos amigos que meus itinerários me trouxeram e que hoje são parte fundamental de meus alicerces, meu amigo e confidente Raphael Silveiras que me mostrou que a paz só é possível através da contradição e caos interior. À Rafaella Garros, minha irmã, que há anos compartilha minhas aflições e medos me auxiliando em toda minha jornada. Também agradeço à minha amiga Natália Novelli que teve participação essencial na minha vida, em uns momentos mais difíceis e maravilhosos, que me mostrou que ainda há esperança para a humanidade, e que há pessoas que são irrevogavelmente boas. Aos amigos e companheiros de trabalho que tive o prazer de conviver ao longo dos meus 10 anos de trabalho no laboratório de Evolução Molecular, Elaine Silva Dias, Adriana Granzotto, Wellington Silva, Elias Carnelossi, Raduan Soleman, Marjorie Silva, Luis Gustavo Galego, Marcelo Jurado, Guilherme Matheus, Edoardo Estevam Lobl, Izabella Luisa Tambones, Lucas Moreira, Bianca Manfré, Felipe Santa-Rosa do Amaral e Camila Vieira, destaco em especial Cecilia Ártico Banho por compartilhar ideais e me apoiar durante todo o processo, não me permitindo desistir, também destaco a relação construída com meu amigo e companheiro de luta e proposito de vida William Vilas Boas Nunes, com o qual não compartilho apenas Ciência mas o desejo de um mundo melhor, e por fim, meus principais agradecimentos é para Daniel de Oliveira, amigo que tive o prazer de fazer em Lyon, sem o qual não teria terminado ou desenvolvido essa tese, sem o qual não teria continuado, agradeço a cada conversa e a todas as utopias compartilhas. Daniel, obrigada! A todos da Equipe do TREEP e do LBBE, Matthieu Boulesteix, Marie Fablet, Annabelle Haudry. À Nelly Burlet e Sonia Martinez pela paciência que tiveram comigo, e em especial agradeço as pessoas que tornaram minha estadia na França a mais agradável possível, apesar do momento em que estávamos, Pierre Marin, Angelo Jacquet, Valentina Rodrigues, Camille Mayeux, Alexis e Vicent Mérel. A todos os técnicos, funcionários, e alunos que de alguma forma contribuíram para a finalização deste processo, em especial aos servidores do Pós-Graduação em Biociências. À Universidade Estadual Paulista “Júlio de Mesquita Filho” – Campus de São José do Rio Preto, que participou da minha rotina e da minha vida nos últimos 10 anos, e que contribuiu para minha formação como cidadã e pessoa. À Université Claude Bernard Lyon 1 e ao Laboratoire de Biometrie et Biologie Evolutive, que fizeram parte da maior e mais desafiadora experiência da minha vida, que além da contribuição acadêmica foram parte fundamental na minha formação como ser histórico- social. Agradeço a FAPESP pela concessão auxílio fornecido ao meu grupo de pesquisa, sob o processo 2016/19271-2, Fundação de Amparo à Pesquisa do Estado de São Paulo (FAPESP). Também agradeço a ANR pela concessão auxílio fornecido ao meu grupo de pesquisa, sob o processo 14-CE19-0016, Agence Nationale de la Recherche. O presente trabalho foi realizado com apoio da Coordenação de Aperfeiçoamento de Pessoal de Nível Superior - Brasil (CAPES) - Código de Financiamento 001. RESUMO Elementos de transposição (TEs) são sequências repetitivas de DNA capazes de se mobilizar dentro e entre genomas que devido a essa habilidade, além de efeitos deletérios, podem trazer vantagens ao genoma hospedeiro, tais como o aumento do repertório proteico, ou influenciarem a regulação gênica. O material de estudo desta tese foi o genoma de Drosophila mojavensis, uma espécie cactofílica, nativa de regiões áridas e semiáridas, que utiliza como sítio de alimentação e reprodução um recurso menos rico que os utilizados pela maioria das espécies de Drosophila e, consequentemente, tal hábito deve estar relacionado a adaptações específicas. Como os TEs são fontes importantes de variação genética, particularmente em resposta ao estresse ambiental, como também podem fornecer matéria-prima para o surgimento de genes codificadores de proteínas (PCG) e de RNAs não codificantes (ncRNA), os quais podem assumir funções celulares importantes, um dos objetivos desta tese foi investigar a associação entre TEs e esses genes em espécies cactofílicas, como estratégia para ampliar nossa compreensão sobre a contribuição dos TEs para a evolução adaptativa. Além disso, mesmo sendo a maioria dos genes comuns aos genomas de machos e fêmeas, genes e TEs podem ser expressos diferencialmente entre testículos e ovários. A análise dos transcriptomas das gônadas torna possível investigar como ou se o viés de expressão de genes e de TEs afeta a reprodução de híbridos interespecíficos. A fim de entender melhor essas questões, foi estudado nesta tese o viés de expressão das gônadas de D. mojavensis e D. arizonae e de seus híbridos recíprocos. Neste estudo mostramos que o conteúdo de TEs representa 13,27% do genoma de D. mojavensis e que a proporção de genes que abrigam sequências de TEs dentro dos éxons é de 5,9%, e que a proporção de ordens de TEs dentro destes éxons difere daquela do genoma geral (LINEs: 78,8% > LTR: 18,5% > TIRs: 2,4% > Helitrons: 0,3%). Entre todas as inserções de TEs encontrados nos genes PCG e lncRNA, apenas 19 sequências apresentaram domínios proteicos preservados. Foram identificados vários genes candidatos à exaptação de TEs no genoma D. mojavensis e a maioria deles abrigam sítios de ligação ao DNA, particularmente domínios zinc-finger, destacando o papel dos TEs como fonte de novidades funcionais aos genomas, evidenciando que esta associação entre sítios específicos de ligação de DNA e TEs é recorrente e talvez uma das fontes principais de novidades genômicas. A maioria dos genes e TEs com viés de expressão estão presentes em ambas as espécies parentais, sendo que em testículos a expressão de genes e TEs é cerca de três vezes maior do que em ovários, destacando que a expressão entre as gônadas em Drosophila é conservada e característica do grupo. Contudo, existem diferenças específicas de genes com viés de expressão entre as espécies parentais e entre os híbridos, que podem ter um papel importante na reprodução. Além disso, as famílias de TEs superexpressas em testículos são mais diversas que os genes superexpressos, esse fato pode estar diretamente relacionado com a influência dos piRNAs na regulação, principalmente devido ao fato dos piRNAs presente nos híbridos ter um padrão diferente dos encontrados nas espécies parentais. Palavras-Chaves: Exaptação. Domínios Zf-BED. Genes De LncRNAs. Sex-biased Genes. Expressão Gênica Gonadal. Genes Diferencialmente Expressos. Isolamento Reprodutivo. Drosophila mojavensis. Drosophila arizonae. ABSTRACT Transposable elements (TEs) are repetitive DNA sequences able to mobilize within and between genomes that due to this ability, besides deleterious effects, can bring advantages to the host genome, such as increasing the protein repertoire or influencing gene regulation. The study material of this thesis was the genome of Drosophila mojavensis, a cactophilic species, native to arid and semi-arid regions, which uses as feeding and breeding site that is less rich in resource than those used by most Drosophila species and, consequently, such habits must be related to specific adaptations. As TEs are important sources of genetic variation, particularly in response to environmental stress, and can also provide raw material for the emergence of protein-coding genes (PCG) and non-coding RNAs (ncRNA), which can assume important cellular functions, one of the objectives of this thesis was to investigate the association between TEs and these genes in cactophilic species, as a strategy to broaden our understanding of TEs contribution to adaptive evolution. Furthermore, even though most genes are common to male and female genomes, genes and TEs may be differentially expressed between testes and ovaries. Analysis of the gonad transcriptomes makes it possible to investigate how or if the expression bias of genes and TEs affects the reproduction of interspecific hybrids. In order to better understand these questions, the expression bias of the gonads of D. mojavensis and D. arizonae and their reciprocal hybrids was studied in this thesis. In this study, we show that TE content represents 13.27% of the D. mojavensis genome and that the proportion of genes harboring TE sequences within exons is 5.9%, and that the proportion of TE orders within these exons differs from that of the overall genome (LINEs: 78.8% > LTR: 18.5% > TIRs: 2.4% > Helitrons: 0.3%). Among all TE insertions found in PCG and lncRNA genes, only 19 sequences showed preserved protein domains. Several candidate genes for TEs exaptation were identified in the D. mojavensis genome and most of them harbor DNA-binding sites, particularly zinc-finger domains, highlighting the role of TEs as a source of functional novelties to genomes, evidencing that this association between specific DNA-binding sites and TEs is recurrent and perhaps one of the major sources of genomic novelties. Most sex-biased genes and TEs are present in both parental species, and there are three times more testis-biased genes and TEs than ovary-biased, highlighting that expression between gonads in Drosophila is conserved and characteristic of the group. However, there are specific differences in sex-bias expression between parental species and between hybrids, which may play an important role in reproduction. Furthermore, testis-biased TE families are more diverse than testes- biased genes, this fact may be directly related to the influence of piRNAs on the regulation, mainly due to the fact that piRNAs from the hybrids have a different pattern than those found in the parental species. Keywords: Exaptation. Zf-BED Domains. LncRNA Genes. Sex-biased Genes. Gonadal Gene Expression. Differentially Expressed Genes. Reproductive Isolation. Drosophila mojavensis. Drosophila arizonae. RÉSUMÉ Les éléments transposables (TEs) sont les séquences d'ADN répétitives capables de se mobiliser à l'intérieur et entre les génomes. Grâce à cette capacité, en plus des effets délétères, ils peuvent apporter des avantages au génome de l'hôte, tels que l'augmentation du répertoire de protéines ou influencer la régulation des gènes. Le matériel d'étude de cette thèse était le génome de Drosophila mojavensis, une espèce cactophile, originaire des régions arides et semi-arides, qui utilise comme site d'alimentation et de reproduction une niche moins riche en ressources que ceux utilisés par la plupart des espèces de drosophiles. Par conséquent, ce comportement doit être liées à des adaptations spécifiques. Comme les TE sont des sources importantes de variation génétique, notamment en réponse au stress environnemental, ils peuvent également fournir de la matière première pour l'émergence de gènes codant pour des protéines (PCG) et d'ARN non codants (ncARN), l'un des objectifs de cette thèse était étudier l'association entre les ET et ces gènes chez les espèces cactophiles, comme stratégie pour élargir notre compréhension de la contribution des ET à l'évolution adaptative. En outre, même si la plupart des gènes sont communs aux génomes mâles et femelle, les gènes et les TE peuvent être exprimés de manière différente entre les sexes et en particulier entre les tissus germinaux. L'analyse des transcriptomes des gonades permet d'étudier comment ou si le biais d'expression des gènes et de TEs affecte la reproduction des hybrides interspécifiques. Afin de mieux comprendre ces questions, le biais d'expression des gonades de D. mojavensis et D. arizonae et de leurs hybrides réciproques a été étudié dans cette thèse. Dans cette étude, nous montrons que le contenu en TE représente 13,27% du génome de D. mojavensis et que la proportion de gènes hébergeant des séquences TE au sein d'exons est de 5,9%, et que la proportion des ordres de TE au sein de ces exons diffère de celle du génome global (LINEs: 78,8% > LTR: 18,5% > TIRs: 2,4% > Helitrons: 0,3%). Parmi toutes les insertions de TE trouvées dans les gènes PCG et lncRNA, seules 19 séquences ont montré des domaines protéiques conservés. Plusieurs gènes candidats pour l'exaptation de TEs ont été identifiés dans le génome de D. mojavensis et la plupart d'entre eux contiennent des sites de liaison à l'ADN, en particulier des domaines à doigts de zinc, soulignant le rôle des TEs comme source de nouveautés fonctionnelles pour les génomes, prouvant que cette association entre des sites spécifiques de liaison à l'ADN et des TEs est récurrente et peut-être l'une des principales sources de nouveautés génomiques. La plupart des gènes et TEs ont une expression sexe biaisé, avec une expression trois fois supérieures de gènes et TE dans les testicules. Cependant, il existe des différences spécifiques dans l'expression des genes et TE, entre les espèces parentales et entre les hybrides, qui peuvent jouer un rôle important dans la reproduction. En outre, les familles de TEs surexprimés dans testicules sont très diversifiées, ce fait peut être directement lié à la production de piRNAs qui sont responsables de la régulation des TE. Mots clés: Exaptation. Domaines Zf-BED. Gènes LncRNAs. Gènes Liées Au Sexe. Expression Gonadique. Gènes Différentiellement Exprimés. Isolement Reproductif. Drosophila mojavensis. Drosophila arizonae. LISTA DE ILUSTRAÇÕES Introdução Figura 1 - Mecanismos pelos quais os TEs podem gerar novidades genômicas ..........29 Figura 2 - Distribuição das espécies de D. mojavensis e D. arizonae utilizadas neste estudo ............................................................................................................................35 Chapter One Figura 1 - Representation of D. mojavensis TE content and percentage of TEs contribution to the mobilome …………………………………………………………………44 Figura 2 - Distribution of lengths of TE insertions into 5´UTR (green), 3´UTR (red), exons (blue), and introns (gray) ……………………………………………………………………..46 Figura 3 - Proportion of TE insertions of each order in the D. mojavensis genes and genome. 5’ UTR (red), 3’ UTR (green), exon (blue), intron (grey) and genome (black) 47 Figura 4 - Hopseu2 insertions within genes of D. mojavensis …………………………...51 Figura 5 - Galileo insertions in the 5’ UTR of multidrug resistance-associated protein 7 gene and the coverage of expression in D. mojavensis transcriptome …………………52 Figura 6 - Homo 10 insertion in exon of lncRNA genes and the coverage of expression in D. mojavensis transcriptome ……………………………………………………………...56 Figura 7 - Additional file 3. Amino acid alignment of exon 1 of the metalloprotease …..76 genes …………………………………………………………………………………………...78 Figura 8 - Additional file 4. Phylogenetic analysis (maximum likelihood) of exon1 of genes that present ZnMc domain …………………………………………………………...79 Figura 9 - Additional file 5. The alignment between the Hopseu2 sequence described and identified in D. pseudoobscura genome (NC_046679.1) and Hopseu2 present in intergenic region of D. mojavensis genome ………………………………………………..80 Chapter Two Figura 1 - Box plot. Box Plot representation of the distribution of log2 fold-change of the parental species and their hybrids ………………………………………………………….91 Figura 2 - Scatter plots of RNA-seq normalized counts of genes of gonads in D. mojavensis and D. arizonae …………………………………………………………………92 Figura 3 - Venn diagram of testis-biased and ovary-biased orthologous genes in D. mojavensis and D. arizonae and reciprocal hybrids ………………………………………93 Figure 4 - Gene enrichment analysis for testis-biased genes in D. mojavensis and D. arizonae and reciprocal hybrids ……………………………………………………………..96 Figura 5 - Ovary-biased genes in D. mojavensis and D. arizonae and reciprocal hybrids. The plot represents the Gene enrichment analysis of all species ……………………….99 Figura 6 - Scatter plots of log10 TE normalized counts in gonads of D. mojavensis and D. arizonae and their reciprocal hybrids ………………………………………………….101 Figura 7 - Venn diagram of testis-biased TEs in D. mojavensis and D. arizonae and reciprocal hybrids ……………………………………………………………………………102 Figura 8 - Relationship between piRNA normalized counts between gonads in D. mojavensis/D. arizonae and their hybrids ………………………………………………...105 Figura 9 - Supp Figure 1 Heatmap. Heatmap of normalized counts of the parental species and their hybrids ……………………………………………………………….…..120 Figura 10 - Supp Figure 2. Difference between gonad replicates ……………….……..121 Figura 11 - Supp figure 3. Principal component analysis of normalized counts of transcriptome gonads replicates of D. melanogaster ……………………………………122 LISTA DE TABELAS Chapter One Tabela 1 - TE insertions in exons, 5´and 3´UTRs of PCGs and in lncRNA genes. TE families identified, genomic position of TE insertion, domain identified in the TE, predicted domain, and predicted gene function ………………………………………48 Tabela 2 - Additional file 1. The NCBI access of each sequence species used to phylogenetic analyses ……………………………………………………………………..…76 Tabela 3 - Additional file 2. TE content and superfamilies percentage identified in D. mojavensis genome (193.83 Mb) ……………………………………………………………77 Chapter Two Tabela 1 - Total number of sex-biased genes for each sample, number of orthologous from D. melanogaster in parentheses and the mean of log2LC …………………………90 Tabela 2 - GO terms of testis-biased exclusive of one parental species ……………….94 Tabela 3 - TOP15 enriched GO terms of DEG found in just one of the parental species …………………………………………………………………………………………………...98 Tabela 4 - Relation between normalized counts from piRNAs and normalized counts from TEs in each gonad for all comparisons ………………………………………….…..103 Tabela 5 - Relation between normalized counts from piRNAs of ovaries and testes for all comparisons ……………………………………………………………………………....….104 Tabela 6 - Supp table 1. GO terms found for genes male-biased between D. mojavensis gonads ……………………………………………………………………………….………..122 Tabela 7 - Supp table 2. GO terms found for genes male-biased between D. arizonae gonads ……………………………………………………………………………….………..122 Tabela 8 - Supp table 3. GO terms found for genes male-biased between H♀mwri♂ari gonads ……………………………………………………………………………….………..122 Tabela 9 - Supp table 4. GO terms found for genes male-biased between H♀ari♂mwri gonads …………………………………………………………………………….…………..122 Tabela 10 - Supp table 5. GO terms found for genes female-biased between D. mojavensis gonads …………………………………………………….…………………….122 Tabela 11 - Supp table 6. GO terms found for genes female-biased between D. arizonae gonads …………………………………………………………………………….…………..122 Tabela 12 - Supp table 7. GO terms found for genes female-biased between H♀mwri♂ari gonads …………………………………………………………………………122 Tabela 13 - Supp table 8. GO terms found for genes female-biased between H♀ari♂mwri gonads …………………………………………………………………………122 LISTA DE ABREVIAÇÕES DEG – Genes diferencialmente expressos DETE – Elementos de transposição diferencialmente expressos Drosophila arizonae – D. ari Drosophila mojavensis – D. moj GO – Gene Ontology H♀ari♂mwri - D. arizonae ♀ x D. m. wrigleyi ♂ H♀mwri♂ari - D. arizonae ♂ x D. m. wrigleyi ♀ LINEs - Long Interspersed Nuclear Elements lncRNAs - RNAs não codificantes longos ncRNAs - RNAs não codificantes PCA - Principal Component Analysis PCG - Genes codificadores de proteínas piRNA - Piwi-interacting RNA RNAm - RNA mensageiro SINEs - Short Interspersed Nuclear Elements TEs - Elementos de transposição TIRs - Terminal inverted repeat UTRs - untranslated regions zf - zinc finger ZnMc - Zinc-dependent metalloprotease SUMÁRIO 1 INTRODUÇÃO ...………………………...................................…………………...........23 1.1 Elementos de transposição e o papel de dessas sequências na evolução genômica .............................................................................................................. .. 24 1.2 Expressão diferencial de genes e TEs em gônadas............................................29 1.3 D. mojavensis e D. arizonae como modelos de estudo.......................................34 2 OBJETIVOS ......................................................................................................... ....37 2.1 Objetivo Geral .........................................................................................................37 2.2 Objetivos Específicos ............................................................................................37 3 Chapter One Transposable elements in Drosophila as drivers of genome evolution: the mobilome of Drosophila mojavensis …………………………………….……….…39 4 Chapter Two Sex-biased genes and transposable elements expression in gonads of Drosophila mojavensis and D. arizonae and their reciprocal hybrids………....…..…81 5 DISCUSSÃO ...…………………..………………………………………………..…….123 6 CONCLUSÕES ...…………………………………………………………………….….131 23 1 INTRODUÇÃO Ao longo das várias décadas de estudos genéticos e genômicos muitas hipóteses foram criadas, validadas ou refutadas, dentre elas, a ideia de que os genes eram unidades fixas nos cromossomos e de que o genoma era estável. O conceito de estabilidade genômica foi desafiado no final da década de 40 por Barbara McClintock (MCCLINTOCK 1950) com a descrição de sequências genéticas em milho que podiam mudar de posição nos cromossomos e, dependendo do local onde se inseriam, podiam alterar reversivelmente o fenótipo do grão. A descrição dessas sequências de DNA móveis, e medianamente repetidas, em praticamente todos os organismos, denominadas elementos de transposição (TEs), iniciou um novo capítulo e um novo universo dentro do corpo de conhecimento sobre a estrutura e funcionamento do genoma. Entretanto, o reconhecimento sobre o papel desempenhado pelos TEs e o impacto de sua existência nos genomas sempre foi uma pauta muito controversa (BIÉMONT, 2010). Atualmente, há consenso que devido à capacidade intrínseca de se movimentarem no genoma grande parte das inserções de TEs tendem a ser desvantajosas, ou neutras; dependendo do local de inserção e da restrição funcional do mesmo, entretanto, eventualmente, algumas inserções podem desempenhar um papel funcional e adaptativo (BROSIUS, 1991; MCDONALD, 1993; MCDONALD, 1995; SHAPIRO, 1999; KIDWELL; LISCH, 2001). Isso é possível porque os TEs são sequências complexas, portadoras de domínios proteicos e sítios de fatores de ligação, que podem ser recrutadas e utilizadas pelo genoma hospedeiro (WICKER et al., 2007; SUNDARAM et al., 2014). Esta fração genômica repetitiva é coletivamente referida como mobiloma (do Inglês, mobilome), juntamente com outras repetições genômicas (KOONIN; WOLF 2008). 24 1.1 Elementos de transposição e o papel de dessas sequências na evolução genômica Devido à grande diversidade de sequência dentro e entre genomas dos diferentes organismos, os TEs foram classificados em diferentes níveis, inicialmente por Finnegan (1989), e posteriormente por Jurka et al. (2005) e Wicker et al. (2007), mas até os dias atuais não há consenso para um sistema universal de classificação (PIEGU et al., 2015). Finnegan (1989) organizou inicialmente a diversidade dos TEs em duas Classes, devido a ocorrência ou não da etapa de transcrição reversa durante o processo de transposição de suas sequências. Foram incluídos na Classe I, os TEs que durante a mobilização apresentam uma etapa de retrotransposição, mediada pela enzima transcriptase reversa que converte a molécula de RNA mensageiro (RNAm) em DNA. Esses TEs, também chamados de retrotransposons, realizam a transposição via um mecanismo denominado copy-and-paste, o qual é caracterizado como transposição replicativa, pois a cada ciclo completo de replicação é produzida uma nova cópia. Por outro lado, os TEs da Classe II utilizam tanto o mecanismo de transposição replicativa (elementos Rolling-Circle), como o mecanismo de transposição sem a etapa intermediária na forma de RNAm. Tal processo de transposição é caracterizado pela excisão do elemento de seu sítio genômico e reinserção em outra região, por um mecanismo denominado "cut-and-paste" (elementos TIRs). Atualizado e ampliado por Wicker et al. (2007), os diferentes sistemas de classificação foram unificados, mantendo-se a divisão de Finnegan (1989) em duas classes e incluindo-se um critério hierárquico que organiza os TEs em classes, subclasses, ordens, superfamílias, famílias e subfamílias. Na Classe I, os retrotransposons foram agrupados 25 em cinco ordens, cujos representantes se diferenciam em dois grupos quanto à presença ou não de longas repetições terminais diretas, as LTRs (do Inglês: Long Terminal Repeats). Três ordens incluem os TEs com LTRs (Retrotransposons com LTRs, elementos DIRS-like e Penelope-like), e duas ordens incluem os TEs sem essas repetições, os non-LTRs (LINEs, do Inglês: Long Interspersed Nuclear Elements, e SINEs, do Inglês: Short Interspersed Nuclear Elements). Os TEs da Classe II se agrupam em duas subclasses: subclasse I (Tc1/Mariner, hAT, Mutator, Merlin, Transib, P, PiggyBac, PIF/Harbinger, CACTA, Crypton) e subclasse II (Helitron e Maverick). Em ambas as Classes existem elementos autônomos, os quais apresentam toda a maquinaria enzimática necessária para realizar sua própria transposição, e também elementos não-autônomos que utilizam a maquinaria dos elementos autônomos para sua movimentação no genoma. Apesar da onipresença de TEs em Eucarya, a abundância dos mesmos é extremamente variável entre os táxons, como em plantas que apresentam uma grande proporção de TEs no genoma, por exemplo em cereais e Zea mays que podem apresentar mais de 85% do genoma constituídos por TEs (LI et al., 2004; SCHNABLE et al., 2009), e o Saccharomyces cerevisiae, cuja porção genômica de TEs corresponde a menos que 4% (CARR et al., 2012). A proporção de TEs também varia entre os grupos de espécies, como em Drosophila, gênero em que o conteúdo genômico de TEs varia entre 4,65% (D. busckii) a 30,80% (D. suzukii) (SESSEGOLO; BURLET; HAUDRY 2016). A maioria das inserções de TEs nos genomas são neutras, ou deletérias, nesse último caso porque as inserções podem interromper a fase de leitura de um gene ou ainda alterar sua transcrição. As inserções deletérias tendem a ser eliminadas dos genomas 26 por seleção purificadora e inserções neutras podem ser eliminadas ou fixadas por deriva genética. Contudo, ocasionalmente, algumas inserções em regiões codificadoras ou em regiões reguladoras tornam-se benéficas, sendo assim, são selecionadas positivamente e mantidas no genoma hospedeiro (FESCHOTTE; PRITHAM 2007). Atualmente, sabe- se que algumas inserções de TEs podem ter um papel adaptativo ao genoma hospedeiro afetando diretamente a evolução genômica e, consequentemente, a evolução das espécies pelo fornecimento de sequências novas tanto codificantes quanto regulatórias (KIDWELL; LISCH 2001, WARREN et al., 2015, LEE et al., 2015, GARCIA-PEREZ et al., CHUONG et al., 2016, CHUONG et al., 2017, JANGAM et al., 2017). Além disso, existem vários processos que podem originar inserções de TEs benéficas, como por exemplo, exonização (Figura 1A) quando uma sequência de TE se insere em um íntron de um gene codificante de proteína gerando novo éxon, uma inserção de TE pode originar um gene codificante de proteína (Figura 1B), ou ainda, atuar como promotor (Figura 1C) ou acentuador (enhancer) (Figura 1D) e, por fim, uma inserção de TE em éxon de genes codificante de proteína (Figura 1E) pode ter parte da sequência recrutada pelo gene. Esse processo de recrutamento de TEs e utilização de tais sequências em novas funções pelo genoma hospedeiro é chamado de exaptação ou domesticação molecular (SINZELLE et al., 2009). Alguns exemplos de exaptação dos TEs nos genomas destacaram o papel adaptativo dessas sequências móveis, como por exemplo o caso dos retrotransposons TART, HeT-A e THARE e os telômeros. Esses LINEs são responsáveis pela manutenção dos telômeros que devido ao processo de retrotransposição desempenham o papel de telomerase, enzima ausente em isentos (PARDUE et al., 2003; SHEEN et al., 1994). Os 27 telômeros em Drosophila abrigam uma enorme quantidade de sequências repetitivas in tandem, resultado de diversos eventos de transposição desses retrotransposons nas extremidades dos cromossomos, locais em que esses LINEs são reversamente transcritos gerando ilhas de TEs com as extremidades 3’ orientadas em direção ao centrômero (PARDUE; DEBARYSHE, 2003). Outro exemplo clássico é o sistema de recombinação V(D)J (TONEGAWA, 1983), que é responsável pelo repertório diversificado de anticorpos e receptores de células T, uma marca registrada da imunidade adaptativa. Este mecanismo, resultante de uma recombinação quimérica entre o transposon Transib e a sequência codificadora do gene RAG1, foi um evento chave na evolução da imunidade adaptativa de vertebrados mandibulados (ROBINSON et al., 2017, KAPITONOV; JURKA 2001). Devido ao advento das técnicas de sequenciamento completo de genomas e das análises de bioinformática, os eventos de exaptação estão sendo cada vez mais investigados e reportados, não apenas aqueles envolvendo genes codificantes de proteínas, como também os genes de RNA não-codificantes (ncRNAs), particularmente os genes de RNA não-codificantes longos (lncRNAs), cujas sequências tem pelo menos 200 nucleotídeos. Esses genes estão diretamente associados como a regulação da expressão gênica, tanto transcripcionalmente quanto pós-transcripcionalmente, interagindo com outras sequências não-codificantes, cromatina e proteínas (revisado em FORT et al., 2021). Estudos recentes têm reportado diversos casos de exaptação em genes lncRNAs. Um dos exemplos mais notáveis, descrito em mamíferos, é a exaptação associada à inativação de um dos cromossomos X, processo que equipara a expressão dos genes do 28 cromossomo X entre os dois sexos, por meio do silenciamento aleatório de um dos cromossomos X das fêmeas, processo esse conhecido como compensação de dose (LYON 1961). Atualmente, sabe-se que esse processo de silenciamento deve-se a várias inserções de TEs exaptados no gene Xist (WUTZ et al., 2002, ELISAPHENKO et al., 2008). Os TEs também são fontes de domínios funcionais para os genes de lncRNAs em vertebrados (JOHNSON; GUIGÓ 2014), como evidenciado em humanos 83% dos genes lncRNAs apresentam pelo menos uma inserção de TE (CARLEVARO-FITA et al., 2019). Em Drosophila, esses genes têm uma estreita relação com desenvolvimento, resposta ao estresse, função gonadal e neurológica (PEASE et al., 2013, LO PICCOLO et al., 2017., MURAOKA et al., 2018, MAEDA et al., 2018, LAKHOTIA et al., 2012, VALANNE et al., 2011). 29 Figura 1. Mecanismos pelos quais os TEs podem gerar novidades genômicas. A) exonização: origem de novo éxon a partir da inserção de TE em um íntron de um gene funcional. B) origem de gene novo: parte de um TE pode ocasionalmente formar um novo gene codificante de proteína. C-D) origem de promotores e enhancers: inserções de TEs na região upstream dos genes codificantes de proteínas podem geram tanto promotores como enhancers novos. E) origem de éxon quimérico: utilização de partes dos TEs em éxons de genes codificantes de proteínas. Fonte: produzida pela autora. 1.2 Expressão diferencial de genes e TEs em gônadas É possível observar dimorfismo sexual entre machos e fêmeas na maioria dos organismos, entretanto, os genomas dos dois sexos são praticamente idênticos; assim sendo, atribui-se o dimorfismo sexual à expressão gênica diferencial entre machos e fêmeas, especificamente aos genes diferencialmente expressos entre as gônadas 30 (HANS; PARSCH 2007, PARISI et al., 2004, ZHANG et al., 2007, CONGRAINS et al., 2018). Tal diferença é significativa e bem documentada em diversos grupos taxonômicos. Em D. melanogaster, por exemplo, três quartos dos genes são expressos diferencialmente entre machos e fêmeas (ASSIS et al., 2012). Esse padrão resulta em diferenças comportamentais, fisiológicas e morfológicas, fundamentando a evolução gênica sexual da expressão gênica, a qual pode ser mais acentuada em tecidos específicos como testículos e ovários. Os genes que são superexpressos em uma gônada em relação a outra são conhecidas como genes com viés sexual (do Inglês, sex-biased genes), especificamente, ovary- biased ou testis-biased (WHITTLE; EXTAVOUR, 2019). Sex-biased genes têm sido vastamente explorados nas últimas duas décadas em diferentes organismos (MEIKLEJOHN et al. 2003; RANZ et al. 2003; ZHANG et al. 2004; ELLEGREN; PARSCH, 2007; HAERTY et al., 2007; WHITTLE et al., 2007; SMALL et al., 2009; ASSIS et al., 2012; WHITTLE; JOHANNESSON, 2013; LIPINSKA et al., 2015). Entretanto, existem poucos estudos sobre genes diferencialmente expressos entre as gônadas. Ao longo dos anos os estudos de sex-biased genes em Drosophila têm mostrado que testis-biased ou male-biased genes, normalmente, apresentam maiores taxas evolução do que ovary-biased ou female-biased genes (ZHANG et al., 2007, ASSIS et al., 2012, MEIKLEJOHN et al., 2002, WHITTLE; EXTAVOUR, 2019). Por outro lado, esse viés pode ser variável entre as espécies, por exemplo, em D. melanogaster a evolução devido a divergência de expressão é maior em male-biased genes, enquanto em D. pseudoobscura esse padrão é observado em female-biased genes (ASSIS et al., 2012). 31 Desde Haldane (1922) se sabe que o isolamento reprodutivo evolui mais rapidamente no sexo heterogamético do que no sexo homogamético, particularmente em insetos. Em Drosophila, nas fases iniciais do isolamento reprodutivo os machos são estéreis e as fêmeas ainda são férteis. Os genes que apresentam maior diferencial de expressão em machos evoluem mais rapidamente que os genes diferencialmente expressos em fêmeas (ASSIS et al., 2012; GRATH; PARSCH, 2012). A alta taxa de evolução dos genes relacionados aos tecidos reprodutivos masculinos poderia levar à maior divergência entre os indivíduos, mesmo sem que tenha ocorrido a esterilidade dos mesmos. A evolução diferencial entre os sexos também pode ser identificada em nível cromossômico. O cromossomo X pode apresentar taxa de evolução maior que cromossomos autossômicos, em populações de tamanho reduzido, pois, em decorrência de deriva genética, mutações no cromossomo X, acabam sendo fixadas randomicamente (MARK et al., 2010; CHARLESWORTH et al., 2018). Assim, devido a potenciais diferenças entre tecidos germinativos e autossômicos é importante analisar separadamente tecidos reprodutivos, a fim de elucidar qual a é relação desses com a evolução do isolamento reprodutivo. Assim como os genes, os TEs também podem ser encontrados diferencialmente expressos entre machos e fêmeas. Ademais, TEs podem ter um importante papel no isolamento reprodutivo, pois atuam nas redes regulatórias de genes e podem também produzir rearranjos estruturais do genoma que afetam o valor adaptativo de híbridos (KIDWELL; LISCH, 2001; SERRATO-CAPUCHINA; MAMUTE, 2018). 32 A expressão dos TEs pode ser alterada tanto por fatores externos, como estresse ambiental (temperatura e patógenos), que podem levar ao aumento dos níveis de expressão (CAPY et al., 2000), bem como por fatores internos, como o estresse genômico decorrente da hibridação intra e interespecífica (revisão em GRANZOTTO; CRUZ, 2015). Em Drosophila, Labrador et al. (1994) mostraram o aumento na transposição do elemento Osvaldo em híbridos de D. buzzatti e D. koepferae, o qual encontra-se reprimido nos genomas das populações parentais. Além disso, já foram relatados eventos de explosão de transposição envolvendo os elementos Helena e Galileo, para o mesmo par de espécies (GUERREIRO, 2014; VELA et al., 2014). Nosso grupo de pesquisa, em uma série de estudos envolvendo as espécies irmãs D. mojavensis e D. arizonae, e seus híbridos produzidos em laboratório, demonstrou a ocorrência TEs diferencialmente expressos em gônadas de fêmeas e machos híbridos em relação a seus respectivos parentais (CARNELOSSI et al., 2014, LOPEZ-MAESTRE et al., 2017, BANHO et al., 2021). Os resultados encontrados indicaram superexpressão de TEs específicos, como os elementos I (CARNELOSSI et al., 2014), GTWIN, Copia e Frogger (LOPEZ-MAESTRE et al., (2017), diferencialmente expressos em relação aos parentais também dependente da direção de cruzamento que produziu esses híbridos. Além disso, em trabalho recente de nosso grupo de pesquisa evidenciou-se a ocorrência de diferentes TEs superexpressos entre machos e fêmeas da prole referente ao mesmo cruzamento interespecífico (BANHO et al., 2021) Ainda, Romero-Soriano et al. (2016) mostraram níveis variáveis de expressão para o retrotransposon Helena entre ovários e testículos em híbridos de D. buzzatii e D. kopeferae, com superexpressão em ovários, e repressão em testículos, da prole F1. Portanto, são relatadas grande diferenças 33 na expressão e regulação de elementos em híbridos, além de ser enviesada em relação ao sexo. Em síntese, o impacto dos TEs nos genomas é significativo e a desregulação dos mesmos pode gerar consequências desvantajosas para os hospedeiros, afetando diretamente o valor adaptativo, especialmente da descendência híbrida, como por exemplo na síndrome de disgenesia híbrida, que foi evidenciada a partir de estudos com o elemento P de D. melanogaster (PICARD, 1976). Essa síndrome é caracterizada pela atrofia gonadal, afetando assim, a viabilidade e fertilidade dos organismos, e ocorre em cruzamentos entre populações maternas com o elemento P ativo ausente (linhagem M) e populações paternas com o elemento P ativo presente (linhagem P). Nesse caso, não ocorrerá no zigoto híbrido a regulação dos elementos P ativos herdados da linhagem paterna, pois a linhagem materna não continha elementos P que ativassem a produção dos piRNAs complementares, e o resultado da mobilização do TE será uma prole onde as fêmeas apresentam os fenótipos da síndrome de disgenesia hibrida (KELLEHER, 2016, HILL et al., 2016). Assim, a regulação dos TEs é extremamente importante para o ambiente genômico, e um dos mecanismos de controle e regulação é a via de regulação dos piRNA. Esse processo tem sido extensivamente estudado, especialmente relacionado a gônadas femininas de Drosophila (KHURANA; THEURKAUF, 2010; JULIANO et al., 2011; BANISCH et al., 2012; HANDLER et al., 2013; WEICK; MISKA, 2014; HUANG et al., 2017). 34 1.3 D. mojavensis e D. arizonae como modelos de estudo Devido ao impacto que eventos de hibridização podem causar no genoma, e ocasionalmente favorecer isolamento reprodutivo, híbridos são ótimos modelos para estudos evolutivos. Espécies irmãs capazes de produzir híbridos em laboratório, como D. mojavensis e D. arizonae, são ótimas candidatas para tais estudos. Essas espécies divergiram muito recentemente, a cerca de 1,5 milhão de anos atrás (SANCHEZ- FLORES et al., 2016), e apresentam isolamento reprodutivo pós-zigótico incompleto e assimétrico, de modo que podem ser produzidos híbridos em laboratório, dependendo da linhagem utilizada, nas duas direções de cruzamento. Ademais, as populações de D. mojavensis apresentam suficiente diferenciação genética, comportamental, ecológica e morfológica que levaram Pfeiler et al. (2009) a classificarem suas populações geográficas em quatro subespécies. Uma dessas subespécies, material de estudo desta tese, D. mojavensis wrigleyi, oriunda de uma população encontrada exclusivamente na ilha de Santa Catalina, na costa da Califórnia (EUA). Também utilizamos como material de estudo nesta tese D. arizonae, cujas populações apresentam ampla distribuição geográfica, mas não são divergentes o suficiente para serem diferenciadas em subespécies (Figura 2). 35 Figura 2. Distribuição das espécies de D. mojavensis e D. arizonae utilizadas neste estudo. Fonte: modificado de BANHO et al., 2021. O genoma de D. mojavensis (GCF_000005175.2), utilizado neste estudo, apresenta 14.714 genes anotados (193.83 Mb). Mais recentemente Petrov et al. (2021), re-sequenciou o genoma de D. mojavensis por meio da tecnologia de nanopore reportando que D. mojavensis possui 15.130 genes (163.17 Mb). Durante a anotação dos 12 genomas de Drosophila (CLARK et al., 2007) reportou-se que o genoma de D. mojavensis apresenta 8,92% de conteúdo de TEs, distribuídos em LTRs (45%), LINEs (28%), DNA transposon (13%) e elementos não identificados (13%). Outros autores também anotaram os TEs do genoma de D. mojavensis. Utilizando como referência o banco de dados de Insecta do Repbase Sessegolo et al. (2016), relataram que o genoma de D. mojavensis abriga ~10% de sequências de TEs e RIUS et al. (2016), utilizando 36 também sequências de novo de uma espécie de Drosophila mais relacionada, D. buzzatii, estimaram quase o dobro desse valor (~ 15%). As espécies irmãs D. mojavensis e D. arizonae pertencem ao grupo repleta de Drosophila, e a maioria das espécies desse grupo utilizam tecidos de cactos em fermentação como sítios de alimentação e oviposição (RUIZ; HEED, 1988a; RUIZ; HEED; WASSERMAN, 1990). A utilização de cactos por espécies do grupo repleta representa uma transposição de desafios ecológicos, uma vez que os cactos têm menor teor de nitrogênio e fósforo do que os outros recursos utilizados pela maioria das espécies de Drosophila, consequentemente, tal hábito deve ser resultado de adaptações específicas (ETGES; ETGES, 1989; MARKOW et al., 2000). Ademais, em cactos, existe uma gama de substâncias alquímicas, tóxicas para a maioria das espécies de insetos. Além disso, os cactos ocorrem em regiões áridas e semiáridas, que impõem dessecação e estresse térmico aos organismos a eles associados (GUILLÉN et al., 2015; MATZKIN; MARKOW, 2009). Assim sendo, as espécies de Drosophila que usam esse recurso passaram por vários desafios ecológicos. Assim, a adaptação das espécies cactofílicas ao uso de hospedeiros com diferentes características bioquímicas e a ambiente xéricos, pode estar relacionada a ocorrência de eventos de exaptação de sequências de TEs que podem ter aumentado a variabilidade do repertório transcricional e proteico. Além disso, como D. mojavensis e D. arizonae produzem híbridos em laboratório, a utilização desses genomas híbridos permite investigar mudanças na expressão gênica e de TEs associadas ao estresse genômico decorrente da hibridização, como também investigar o papel da expressão diferencial entre as gônadas no isolamento reprodutivo pós-zigótico. 37 2 OBJETIVOS 2.1 Objetivo Geral Este estudo teve três eixos centrais e complementares. O primeiro foi caracterizar o mobiloma de D. mojavensis, visando a compreensão da abundância e diversidade do conteúdo de TEs do genoma para identificar possíveis eventos de exaptação. O segundo eixo foi avaliar a extensão do viés de expressão de genes e TEs em gônadas das espécies irmãs D. mojavensis e D. arizonae. O terceiro eixo foi avaliar se há viés de expressão de genes e TEs entre as gônadas dos híbridos, e se essas sequências estão relacionadas a funções reprodutivas, sendo assim são candidatos a genes associados ao isolamento pós-zigótico. 2.2 Objetivos específicos a) Anotar e classificar até o nível de família os TEs de D. mojavensis (GCF_000005175.2) e D. arizonae (GCA_001654025.1); b) Caracterizar a diversidade e abundância do mobiloma de D. mojavensis; c) Identificar possíveis eventos de exaptação em genes codificadores de proteínas e de lncRNAs no genoma de D. mojavensis; d) Realizar a análises de expressão diferencial entre as gônadas de D. mojavensis e D. arizonae e seus híbridos recíprocos e avaliar a proporções de genes e TEs com viés de expressão em ovários e testículos; 38 e) Investigar se os genes e TEs com viés de expressão nas duas gônadas das espécies parentais são diferentes dos encontrados nas gônadas dos híbridos e se o viés está associado ao isolamento pós-zigótico incipiente. f) Inferir se há genes testis-biased associados a funções reprodutivas e se podem ser associados ao isolamento pós-zigótico incipiente; g) Quantificar o conteúdo de piRNAs em ovários e testículos e avaliar sua influência na regulação de TEs nas duas gônadas. 39 3 Chapter One - Transposable elements in Drosophila as drivers of genome evolution: the mobilome of Drosophila mojavensis Note: Submitted to Genome Biology | Print: 1474-760X | Electronic: 1465-6906 Simão MC1,2, Oliveira DS1, Vieira C2 and Carareto CMA1 1UNESP - São Paulo State University, Department of Biology, São José do Rio Preto, São Paulo State (SP), Brazil 2 Université de Lyon, Université Lyon 1, CNRS, Laboratoire de Biométrie et Biologie Evolutive UMR 5558, F-69622 Villeurbanne, France Correspondance: simaomc.bio@gmail.com, claudia.carareto@unesp.br Abstract Background: Transposable elements (TEs) are repetitive DNA sequences able to mobilize within and between genomes, and due to this ability, they can bring advantages to the host as expanding the protein-coding sequence repertoire as well as playing a role in gene regulation. Due to the biological and ecological characteristics of Drosophila mojavensis, which is a cactophilic species capable of surviving in challenging environments, the public genome of this species was chosen to analyze the TE content and the mobilome characteristics aiming to investigate their association with protein-coding genes (PCG) and long non-coding RNA (lncRNA) genes. Results: The TE content of D. mojavensis represent 13.27 % of the genome and this mobilome is composed of 38.9 % of LTR retrotransposons, 22.8 % of TIRs, 21.7% of Helitrons and 16.6 % of LINEs. The most abundant TE families in the genome are Helitrons-1 and Homo-6, which represents 24.7% and 21.7% of total insertions. The proportion of genes that harbor TE sequences within exons was 5.9 % and the proportion of TE orders inside these exons differ from the general genome (LINEs: 78.8 % > LTR: 18.5 % > TIRs: 2.4 % > Helitrons: 0.3 %). From the total amount of TE insertions within PCG and lncRNA genes, only 19 sequences carried preserved protein domains. Conclusions: Our analysis showed 14 genes candidates to TE exaptation in D. mojavensis genome. All the TE sequences inserted into exons of PCG (6) or of lncRNA genes (8) harbor DNA- binding sites, mainly zinc binding sites, evidencing a role of TEs as source of functional domains and showing that this association between specific binding sites and TEs is recurrent and maybe a main source of genomic novelties. Keywords: Exaptation. Zf-BED Domains. LncRNA Genes. 40 BACKGROUND Transposable elements (TEs) are repetitive DNA sequences that have the ability to mobilize within and between genomes and due to this ability challenged the paradigm of a static genome early in the 1950s, as proposed by Barbara McClintock [1]. This repetitive genomic fraction is collectively referred to as mobilome, together with other genomic repeats [2]. The TE content is quite variable in eukaryotes, both quantitatively and qualitatively. In their broader organization TEs are classified hierarchically into two classes, based on their transposition mechanism. Class I TEs (Retrotransposons) have a copy-and-paste mechanism in which the TEs are reverse transcribed and the complementary DNA (cDNA) is inserted elsewhere in the genome [3,4]. According to the hierarchical classification proposed by Wicker et al. [5], the retrotransposons are grouped into five orders, whose representatives are classified into two groups according to the presence of long terminal repeats (LTRs) or their absence (non-LTRs). Three orders include TEs with LTRs (LTR retrotransposons, or just LTRs, DIRS-like and Penelope-like elements), and two include non-LTR retrotransposons (LINEs: Long Interspersed Nuclear Elements and SINEs: Short Interspersed Nuclear Elements). Class II TEs have a transposition mechanism that directly cut the DNA molecule, and can be grouped into two Subclasses (Subclass I: TIR and Crypton and Subclass II: Helitrons and Maverick) that differ by the number of DNA strands cleaved during the transposition process, resulting in either conservative transposition (cut-and-paste), in which both strands are cleaved, and the excised sequence inserts into a new site, or replicative transposition (copy-and-paste), in which only one strand is cleaved, involving its strand displacement and formation of a rolling-circle structure, originating a second TE copy [6]. 41 TEs are virtually present in all eukaryotic genomes analyzed so far [7 – 9]. However, their distribution within and among genomes depends on the organism and it is quite variable among eukaryotic genomes [10]. For example, the genome of Belgica antarctica, an insect endemic in Antarctica, has less than 1% of TEs [11], and the Australian lungfish Neoceratodus forsteri contains 90% of TEs [12]. In insects, which is a very diverse group, large differences in TE content are observed, even within close related species. For example, in mosquitoes, Aedes aegypti has 55% of its genome composed by TEs, while in Anopheles gambiae the genomic content is about ~20% [13]. In Drosophila, the TE content is also quite variable, as in some species in which the TE content is low, as in D. busckii (4.65%) and in others it is significantly higher, as in D. suzukii (30.80%) [14]. Drosophila mojavensis is a cactophilic species that belongs to the repleta group (Drosophila subgenus) that has been widely used in speciation studies. However, its mobilome and the potential role of TEs PCG and lncRNA evolution is noticeably understudied. The few reports found in the literature concern the total TE content and description of their classes and orders, as in the first annotation published in the Drosophila 12 Genomes Consortium [15], or later, by Rius et al [16], which annotated 8.9 % and 15.35 % of TEs, respectively, in the D. mojavensis genome. As it has been widely reported over the past 40 decades, that TEs have been considered genomic parasites because of their ability to replicate faster than the host genome and to cause deleterious effects when they insert into coding regions, interrupt open reading frames, or affect gene expression [17 - 20]. However, in recent decades, it has become a consensus that TEs play a fundamental role in the evolution of living beings. 42 Indeed, the insertion of TEs, despite the potential deleterious consequences, can influence gene expression and function, increasing the variability of the transcriptional and protein repertoire [21 - 25]. Also, exaptation of TE sequences can occur when their sequences are co-opted by the host genome [26, 27]. Such changes can generate new genes and functions that have a key role in adaptation and speciation [26, 28, 29]. Thereby, an unexpected TE insertion into a coding or regulatory region can create evolutionary innovations. In the last two decades several exaptation examples have been reported in various organisms. The classical example in Drosophila is the exaptation of three LINEs, HeT-A TART and THARE, which play the role of telomerases, the enzymes responsible for the maintenance of telomeres throughout cell cycles, which is absent in Drosophila and other insects [30, 31]. Another classical example is the V(D)J recombination process [32], a mechanism that results in a diverse repertoire of antibodies and T cell receptors, which is a hallmark of the adaptive immune system appearance. This mechanism is associated with an exaptation of Transib, a DNA transposon, by the RAG1 gene, resulting in a chimeric recombination between the TE and the host gene coding sequences, and was a key event in the evolution of jawed vertebrate adaptive immunity [6, 33]. Exaptation events have also been described in lncRNA genes. These genes regulate gene expression both transcriptionally and pos-transcriptionally and interact with chromatin, proteins and other coding and non-coding RNAs. For example, in mammals, the X chromosome silencing process is associated with the exaptation of several TEs in six exons of the lncRNA gene Xist [34, 35]. In addition to the classical example of Xist, other studies have reported close association of lncRNA domains with transposable 43 elements (reviewed in Fort et al [36]), which lead Johnson and Guigó [37] to propose that TE insertions within lncRNA genes provide a source of functional domains for the transcripts. The Drosophila lncRNA genes are strongly related to organ development, stress response, gonadal function or neurological function [38 – 43] and the association between their domains and TEs is not well known. In summary, the role of TEs to provide evolutionary novelties has accumulated in recent decades. In this study we analyzed the D. mojavensis mobilome, as a model of species exposed to stressful environments, seeking to identify TEs that could be sources of genome innovations. Contrary to most of drosophilids, which mainly use fermented fruits from plants in humid forests, this and other cactophilic species face major ecological challenges as they live in arid and semiarid environments and use cacti as feeding and breeding sites, which in addition to having less nitrogen and phosphorus than the resources used by most other Drosophila species, also produce alchemical substances that are toxic to most insects [44]. So, ecological changes in the use of cacti are believed to have involved a large number of adaptive changes in the physiology, reproductive biology and behavior of cactophilic species [45]. Therefore, due to the increase in the variability of the transcriptional and protein repertoire that can be generated by exaptation events, we analyzed possible genes and TEs that can be candidates for these events in the cactophilic species D. mojavensis. Here, we provide details of the D. mojavensis mobilome and evidenced influence of TEs to genome evolution of this species. We found 14 genes candidates for harboring TEs exapted, six having TE sequences inserted into PCGs and 8 into lncRNA genes, most of them associated to zinc motifs. 44 RESULTS Proportion of transposable elements in the D. mojavensis genome We annotated 46,471 TEs, which correspond to 13.27% (27.03 Mb) of D. mojavensis genome (193.83 Mb). From this content (Figure 1A), the mobilome is constituted of 55.49 % Class I (38.87 % are LTRs and 16.63 % are LINEs) and 44.51 % of Class II elements (22.8 % are TIRs and 21.7 % are Helitron), plus few unknown repeats (Figure 1 B and Additional file 2). Among the superfamilies, 287 families were identified and most abundant are Helitron-1 (11,496 insertions) and Homo-6 (9,857 insertions), respectively from Helitron and the hAT superfamilies, which together represent 45.86 % of insertions in the D. mojavensis mobilome (Helitron-1: 24.69 % and Homo-6: 21.17 %) (Figure 1C). Figure 1. Representation of D. mojavensis TE content and percentage of TEs contribution to the mobilome. (A) TE contribution to the genome, in base pairs; (B) TE orders contribution relative to the total TE fraction in the mobilome; (C) Total TE copy numbers of all TE superfamilies identified in the genome (GCF_000005175.2_ moj_caf1). Source: produced by the author. 45 Association between transposable elements and D. mojavensis genes The release GCF_000005175.2_ moj_caf1 of the D. mojavensis genome has 14,714 annotated genes, of which 13,329 (90.5 %) are PCGs, 979 (6.6 %) are lncRNA genes and 404 (2.7 %) are other ncRNA genes. TEs were annotated in four regions of the D. mojavensis genes: (i) exons, (ii) 5´and 3´ UTR and (iii) introns. From the 46,471 TE insertions, 1,722 sequences (3.8 %) were found within exons, 77 in 5’ UTR (0.16 %), 63 in 3’ UTR (0.13 %) and 7,126 (15.3 %) in introns. There are 872 PCGs (5.9 %) and 36 lncRNA genes (0.24 %) that have at least one TE insertion in exons. We found 52 3’UTRs and 63 5’UTRs harboring at least one TE insertion which correspond to 0.35 % and 0.42 % of D. mojavensis genes, respectively. Finally, 1,528 (10.38 %) genes harbor at least one TE insertion in their introns. The median length of TE sequences within exons is 168 bp, in the 5’ UTR is 125.5 bp, in the 3’ UTRs is 148 bp and in the introns is 265 bp (Figure 2) which shows that most of these insertions are very small. The TE median of insertion lengths was significantly different between the gene regions (Kruskal-Wallis: χ2 = 607.28, df = 1, p < 2.2e-16). 46 Figure 2. Distribution of lengths of TE insertions into 5´UTR (green), 3´UTR (red), exons (blue), and introns (gray). Source: produced by the author. The orders of TE insertions in the exons are highly heterogeneous and do not correspond to the observed proportion in the genome (Figure 3, χ2= = 5082.1, df = 3, p- value < 2.2e-16). Although TIR and Helitron are the most frequent TE orders in the genome, within the exons the LINE order is predominant (78.8%) as well as in the UTRs, ranging from 62.3 % (5’ UTR) to 42.8 % (3’ UTR). 47 Figure 3. Proportion of TE insertions of each order in the D. mojavensis genes and genome. 5’ UTR (red), 3’ UTR (green), exon (blue), intron (grey) and genome (black). Source: produced by the author. Transposable element insertions containing a protein domain in gene regions From the 1,722 TE insertions found within exons, only 14 presented protein domains. From these, six were PCGs and eight were lncRNA genes (Table 1). Among them, 12 DNA transposons and two LINEs (LOA and Jockey) were identified. From the DNA transposons, two superfamilies were identified (hAT superfamily and P superfamily). In the UTR regions two insertions were found in the 5’UTRs (from Mariner and Bel-Pao superfamilies) and three in the 3’UTR (P superfamily). 48 Table 1 TE insertions in exons, 5´and 3´UTRs of PCGs and in lncRNA genes. TE families identified, genomic position of TE insertion, domain identified in the TE, predicted domain, and predicted gene function. Source: produced by the author. TE families Class – Order-Superfamily Genomic position Gene Domain name Predicted domain function Predicted gene function TE length Protein-coding genes exons Hopseu 2 Class 2 - DNA transposon – hAT superfamily NW_001979112.1-31985820- 31986148 LOC26528709 ZnMc superfamily digest the egg envelope or chorion egg envelope or chorion digestion 329 Hopseu 2 Class 2 - DNA transposon – hAT superfamily NW_001979112.1-20457716- 20458496 LOC6574081 ZnMc superfamily digest the egg envelope or chorion egg envelope or chorion digestion 780 Hopseu 2 Class 2 - DNA transposon – hAT superfamily NW_001979112.1-10652216- 10652542 LOC6573082 ZnMc superfamily digest the egg envelope or chorion egg envelope or chorion digestion 326 Hopseu 2 Class 2 - DNA transposon – hAT superfamily NW_001979112.1-10656673.- 10657001 LOC6573084 ZnMc superfamily digest the egg envelope or chorion egg envelope or chorion digestion 328 Loa-1 Class 1 – LINE- Loa superfamily NW_001979112.1-23748677- 23750117 LOC6574428 DUF4780 superfamily uncharacterized uncharacterized 422 Jockey-3 Class 1 - LINE - Jockey NW_001979118.1-3497661- 3498037 LOC6585671 PRK10263 superfamily DNA translocase putative RNA exonuclease pqe-1 377 Protein-coding genes 5´UTR Galileo Class 2 - DNA transposon - P superfamily NW_001979113.1-5521622- 5522362 LOC6576030 THAP superfamily DNA-binding P-element multidrug resistance- associated protein 7 741 Galileo Class 2 - DNA transposon- P superfamily NW_001979113.1-5524088- 5524884 LOC6576030 THAP superfamily DNA-binding P-element multidrug resistance- associated protein 7 797 Protein-coding genes 3´UTR Mariner-12 Class 2 - DNA transposon – Mariner superfamily NW_001979123.1-2251566- 2251991 LOC6586721 HTH_Tnp_Tc3_2 superfamily transposase uncharacterized 426 Mariner-18 Class 2 - DNA transposon – Mariner superfamily NW_001979121.1-914481- 914784 LOC6586201 rve superfamily integrase zinc finger protein 2 homolog 304 BEL-1 Class 1 – LTR - Bel-Pao NW_001979112.1-1126194- 1126452 LOC6572007 RT_like superfamily reverse transcriptase forkhead box protein P1 259 lncRNA genes Galileo Class 2 - DNA transposon- P superfamily NW_001980364.1-3731-4453 LOC26527540 THAP superfamily DNA-binding P-element uncharacterized 723 Galileo Class 2 - DNA transposon- P superfamily NW_001979113.1-5514104- 5515009 LOC26528160 THAP superfamily DNA-binding P-element uncharacterized 802 Homo10 Class 2 - DNA transposon – hAT superfamily NW_001979113.1-31971965- 31972150 LOC26527136 zf-BED superfamily transposase binding uncharacterized 186 Homo10 Class 2 - DNA transposon – hAT superfamily NW_001979119.1-1869453- 1869638 LOC26527196 zf-BED superfamily transposase binding uncharacterized 186 Homo10 Class 2 - DNA transposon – hAT superfamily NW_001979119.1-2610303- 2610451 LOC26527472 zf-BED superfamily transposase binding uncharacterized 149 Homo10 Class 2 - DNA transposon – hAT superfamily NW_001979117.1-535757- 535942 LOC26527778 zf-BED superfamily transposase binding uncharacterized 186 Homo10 Class 2 - DNA transposon – hAT superfamily NW_001979113.1-6896655- 6896840 LOC26528070 zf-BED superfamily transposase binding uncharacterized 186 Homo5 Class 2 - DNA transposon – hAT superfamily NW_001979121.1-294278- 294426 LOC26528413 zf-BED superfamily transposase binding uncharacterized 149 49 Insertions into protein-coding genes DNA transposon - Hopseu2 Hopseu2 is a TIR transposon, from the hAT superfamily, which has five insertions in the D. mojavensis genome. Four of them are within the first exon of the genes coding for metalloproteases, the LOC26528709 and LOC6573082 encode a Low choriolytic enzyme, LOC6574081 encodes a zinc metalloproteinase nas-6 and LOC7573084 encodes an astacin (Fig. 4). These insertions have between 300 – 700 bp and lack transposase domains but present a ZnMc domain (Zinc-dependent metalloprotease), which belongs to a group of zinc-dependent proteolytic enzymes with a HExxH zinc-binding active site. Members of this domain family may have an amino terminal propeptide, which is cleaved to produce the active protease domain [46]. The fifth Hopseu2 insertion found is intergenic and the only that presents a transposase domain. The search made in Drosophila genomes to investigate the range and homology of these metalloprotease genes returned hits in 22 Drosophila genomes (Additional file 3). The ML phylogeny of exon 1 of these genes resulted in a well-supported clade (100 %) that clusters the sequences of species from the Sophophora subgenus that is relatively congruent with the species phylogeny (Additional file 4), which is an internal clade of a bigger one that includes sequences of species from three different subgenera (Drosophila, Lordiphosa and a sequence of D. willistoni, which belongs to the Sophophora subgenus), also highly supported. This paraphyletic clade is also highly supported (98%). There is also in the phylogeny a set of 10 sequences from three species of the repleta group (D. hydei, D. arizonae and D. mojavensis) that forms two internal clusters with high support 50 (71 % and 99 %); in each one, genes from D. mojavensis and D. arizonae genomes can be found. Three groups of metalloproteases in D. mojavensis and D. arizonae can be observed clustering together, suggesting that they are paralog genes, probably due to a recent duplication in the common ancestor of these species. A search for expression of the metalloprotease genes made in an illumina RNAseq transcriptome aggregate available on NCBI showed that the transcription of all these genes is not affected by the TE insertions (Figure 4). However, the expressions levels of these genes are different. These chimeric transcripts of Hopseu2 in the D. mojavensis transcriptome reinforce our hypothesis of a TE exaptation. This hypothesis is also reinforced by the presence of a Hopseu2 sequence carrying the transposase dimerization domain found in the intergenic region of the D. mojavensis genome (scaffold NW_001980292.1, position 4591-5521), which is 64% similar (Additional file 5) to Hopseu2 described in D. pseudoobscura [47], but it has no similarity with the five Hopseu2 sequences inserted into the genes. 51 Figure 4. Hopseu2 insertions within genes of D. mojavensis. The “>” means the sense orientation of the gene. The Hopseu2 insertions and RNA-seq data were obtained in the GCF_000005175.2_ moj_caf1 release of D. mojavensis genome available on NCBI. Source: produced by the author. DNA transposon - Galileo Three insertions of Galileo (ranging from 741 - 797 bp), which belong to the P superfamily and Galileo family, were found in tandem in the 5´UTR of multidrug resistance-associated protein 7 gene, one of them is in the intron of this UTR (Figure 5). All insertions, which are in tandem in the UTR, presents the THAP domain of Galileo, which is involved in site-specific DNA binding function and has a highly conserved zinc finger motif, characteristic of the P superfamily [48]. The multidrug resistance-associated protein 7 gene belong to one subfamily of the ABC transporter family, which is subdivided into seven subfamilies, and some of the proteins they encode are involved in the transport of endogenous and xenobiotic components [49]. The multidrug resistance-associated 52 protein 7 gene have two transcripts and both present at least one Galileo in the sequence. These insertions were not found in the multidrug resistance-associated protein 7 gene in other Drosophila genomes from the repleta group, except in the D. navojoa genome, which belongs to the D. mojavensis complex with D. arizonae (LOC108658602), which present three Galileo insertions in upstream region of these gene. Figure 5. Galileo insertions in the 5’ UTR of multidrug resistance-associated protein 7 gene and the coverage of expression in D. mojavensis transcriptome. The Galileo insertions and RNA-seq data were obtained in the GCF_000005175.2_ moj_caf1 release of D. mojavensis genome available on NCBI. Source: produced by the author. DNA transposon - Mariner Two DNA transposons from the TIR order (Tc1-Mariner superfamily and Mariner family) were found in the 3’ UTR of two genes. The first, with 304 bp, is inserted in the 3’ UTR of the zinc finger gene LOC6586721, presents the integrase domain (rve 53 superfamily) from the transposase, which mediate the DNA integration into the genome. The second insertion, with 426 bp, is within a 3´ UTR of a gene (LOC6586721) with uncharacterized function. This insertion presents the HTH_Tnp_Tc3_2 superfamily domain, which is a catalytic domain characteristic of Mariner and also related to specific DNA binding. Retrotransposon insertions Three retrotransposons are inserted in PCGs, two LINEs from LOA-1 and Jockey- 3 families and the LTR retrotransposon Bel-1. The LOA-1 insertion, with 422 pb, was found within the only exon of an uncharacterized gene (LOC6574428) and it presents the DUF4780 superfamily domain that has an uncharacterized function and is widely found in Eucarya [50]. The wide distribution of this domain suggests that it plays a functional role. The Jockey-3 sequence, with 377 bp, is inserted within the first exon of the gene LOC6585671, which has a function of RNA exonuclease, and the domain within the TE insertion has a DNA translocase function. The gene function and the TE domain function are characterized as basic molecular function for cellular replication. The BEL-1 sequence, with 259 bp, was inserted into the 3’ UTR of the gene LOC6572007, which encodes for a forkhead box protein P1 gene and present the domain RT-like superfamily. 54 DNA transposon insertions in lncRNA genes Galileo Three DNA transposon Galileo were found in lncRNA genes, one of them within the first exon of the gene LOC26527540, with also has one insertion within third intron, and one in the first exon of LOC26528160 gene. Both genes present more than one expressed transcript, as found in the RNA-seq data of D. mojavensis available. The length of all Galileo insertions in these genes are similar ~ 723-756 bp, which is much smaller than a complete sequence found in the D. mojavensis genome (~ 6,576 bp) and all of them present the Galileo characteristic domain THAP (THanatos-Associated Protein) domain intact. Homo Homo are DNA transposons that belong to the hAT superfamily (Homo family). Six Homo insertions were found within lncRNA genes, of which five belong to Homo10 family and one to Homo5. Two insertions have 149 and four have 186 bp and all of them harbor a complete zf-BED domain (zinc finger BED domain) of ~160 bp (between positions 35 and 164 of the insertion). The lncRNA genes, which are ~800 bp long, have three exons and the TE sequences are within the third exon in all of them. Figure 6 represents the structure of these genes using LOC26527136 as example, which has a Homo10 insertion of 186 bp. These genes are highly conserved with a genetic divergence (p-distance) of 55 0.06%. This similarity suggests that these genes result from duplication events. Interestingly, that high similarity of lncRNA genes is restricted to lncRNA genes that have the ZBED domain, and these six genes do not represent the pattern of similarity found among lncRNA genes of D. mojavensis, since the remaining 973 lncRNA genes present non-alignable sequences. Figure 6. Homo 10 insertion in exon of lncRNA genes and the coverage of expression in D. mojavensis transcriptome. The insertion was obtained in the last release of D. mojavensis genome available on NCBI. Source: produced by the author. DISCUSSION Due the TE importance as a source of genomic novelties [51 - 53], several studies have highlighted the role of these mobile sequences in adaptation [54 – 56], particularly in species that challenged biological barriers as well as those that colonized and inhabit adverse environments [57, 58]. D. mojavensis, the subject of this study, is a species native to arid and semiarid environments and specialized in the use of cacti as feeding and breeding sites and that contain compounds toxic to most insects. Because of these 56 peculiarities, the genome of D. mojavensis is interesting to understand better the role of TEs in organismic adaptation. TE content in D. mojavensis mobilome, total and in gene regions To access the TE insertions in the coding regions of the D. mojavensis genome, a preliminary annotation of its mobilome was made and the repeats to the family level were classified, which is an original contribution of this study. The first annotation of the D. mojavensis mobilome was that of the Drosophila 12 Genomes Consortium [15] and since then only two studies made a more complete annotation from this genome at the order or superfamily levels [14]. The TE content presented in our study represents 13.27% of the D. mojavensis genome, very similar to the 14.35% reported by [16] and our classification of the TE is also in agreement to the one reported by this author with 57.3% and 42.6% and 53.6% and 43.1%, respectively of Class I and Class II. These values differ from Sessegolo et al. [14] in which Class II is more abundant than Class I. These annotations resulted in different contributions of the TE orders to the genome content, as for example when comparing our data LTR > TIR > Helitron > LINE and the Sessegollo et al. [14], in which TIR > LTR > LINE > Helitron. This difference occurs also comparing our data with that of Ruis et al. [16] in which Helitron ~ LTR > LINE > TIR. These differences are not unusual in such studies and depends on several factors such as the tools and the TE library used to annotate the repeats, as well as the sequencing technology, which impacts the quality of genome assemblies [16, 59 – 61]. So, as it has been widely acknowledged, a direct comparison of TE annotation may be risky. Even so, it is necessary to point out 57 that our objective was not to make a detailed study of the D. mojavensis mobilome but presenting an overview of the TE insertions up to the family level that would enable the identification of specific insertions into genes, the aim of this study. The documentation of roles of TEs as sources of genetic novelties and their positive impact in genome evolution is growing [25, 31, 58, 62]. It is also well known that the insertion of TEs within PCGs can drastically affect the genomic stability, due to functional constraints of coding sequences; thereby, TE insertions tend to be eliminated from these regions by purifying selection [63 - 65]. The number of TE insertions within these regions is thus expected to be lower than in other genomic regions. In the D. mojavensis genome, 5.9% of genes have at least one TE insertion located at the exons. We found significantly fewer TEs in exons and UTRs than in introns, as expected, but the TE proportion in the exons seems to be high, even if the vast majority of these insertions being very small (median: 168 bp). In D. melanogaster referenced genome, 513 TE insertions are inside exons [65]. However, the comparison with our analysis is not really valuable, since the TE annotated in the reference genome of D. melanogaster only represent 6% of the genome, and we know this content should be up to 12% [14]. For this reason, the results found in D. mojavensis are representing the TEs distribution in a genome-wide manner, possibly overestimated, whereas the data available for D. melanogaster is relying only on bona fide TE insertions. Also, misannotation of either genes or TEs in the D. mojavensis genome is very likely, as this genome is of lower quality than the D. melanogaster one. Even so, most of these insertions were manually checked and confirmed as to their insertion sites, as can be seen with the 18 insertions that kept the protein domains, whose genomic coordinates are shown in Table 1. 58 In order to be most conservative, we have analyzed TE insertion that were detected in exons and performed a more detailed analysis of the ones that presented a catalytic domain. Indeed, most of the well documented examples of exaptation concern DNA transposons [66, 67] probably because of the presence of the N-terminal DNA-binding domain (DBD) and the C-terminal transposase catalytic domain in the transposase sequence. In our data set, more than 95% of TE insertions found within the exons and UTRs were LINEs, followed by LTR retrotransposons, most of which were small, and did not contain protein domains, and only 19 genes have TE insertions with these domains, the great majority were DNA transposons (15) that belonged to only three superfamilies, hAT, Mariner and P. Transposable elements exaptation Exaptation of hAT transposase in eukaryotic organisms have been extensively documented [68 – 71]. Hickman et al. [72] suggested the Zn-binding BED domain located at the N-terminus of the hAT transposase is involved in nuclear localization and non- specific DNA binding, while the proximal C-terminal protein block is related to specific DNA-binding to transposon terminal and sub-terminal repeats. Fusion between domesticated transposases and zinc finger protein domains has also been reported [73]. Moreover, it was proposed that exaptation involving hAT transposons and zinc domains and specific DNA-binding proteins were originated from several and independent hAT exaptation events in different eukaryotic lineages, such as fish and insects [71, 74]. 59 Hopseu2 transposon: Is it a case of exaptation? In this study we identified four copies of Hopseu2, a hAT transposon, inserted into metalloproteases genes of the D. mojavensis genome. This hAT transposon was first identified in the genome of D. pseudoobscura as a single copy with 4,962 bp presenting an intriguing characteristic. It apparently encoded a functional protein of 821 aa, but did not present TIRs (Terminal Inverted Repeats) nor TSDs (Target Site Duplication), suggesting an ancient domestication [47]. The authors argued that ‘‘inactive’’ elements like Hopse2 should show some degeneration in their coding sequence, except if they had been domesticated as a new gene, what could be the case, but demonstration was not given. From the four metalloproteases genes presenting the ZnMc domain, two encode the Low choriolytic enzyme, a protease with catalytic activity described in the inner layer of fish egg envelope [80], as well as associated to egg hatching [75, 76]. The search for the homolog sequences using the exon 1 sequence of Low choriolytic gene, which has a 329 bp fragment identified as Hopseu2, in genomes of four Drosophila subgenus (Sophophora, Drosophila, Lordiphosa and hawaian Drosophila), returned sequences from 22 genomes with high similarity among Hopseu2 insertion and exon 1 of metalloprotease genes in several species of Drosophila. Moreover, the clustering of three genes of D. mojavensis in different clades of the phylogeny, one inside the clade that includes species of four Drosophila subgenus and two clustering in two well- supported clades that includes only species from the mojavensis complex (D. mojavensis, its sister species D. arizonae), show that these genes belong to a gene family containing at least three metalloprotease paralogs. We also showed that these genes with the TE 60 insertion are expressed, suggesting the potential functional role of the Hopseu2. In addition, the observation of a Hopseu2 insertions in intergenic region of D. mojavensis genome, with a transposase domain and with more than 60% of similarity with the Hopseu2 described in D. pseudoobscura, support the suggestion of an exaptation event. Synthesizing, we suggest that Hopseu2 possible exaptation was an ancient event due to three main pieces of evidences: (i) impossibility of transposase identification in the TE sequence within genes, (ii) presence of ZnMc domain in Hopseu2 in the catalytic genes of Drosophila species belonging to four subgenera, and (iii) non-interruption of transcription by the TE insertion in D. mojavensis. We found Hopseu2 domesticated in 17 Drosophila genomes, however, all 27 analyzed Drosophila genomes have ZnMc domain in these catalytic genes. Our hypothesis is that due to the catalytic characteristics of the ancient transposases, these sequences (~100 bp) were recruited by genes as an extension of their own domains. Moreover, this is the first report of one Hopseu2 sequence with a transposase domain in intergenic region and that did not go through the exaptation event, possibly derived from the ancestral sequence. The accumulation of differences between the exapted and genomic sequences reinforces the hypothesis that Hopseu2 sequences are good candidates for exaptation events. 61 lncRNA genes and TEs: a source of possibilities The zf-BED domains - A long lasting history about exaptation ZBED are genes characterized by the zf-BED domain (zinc finger) that belong to a gene family originated from transposase of hAT DNA. This domain was named BED because of Drosophila exaptation event due the recruitment from the DNA transposons in BEAF and DREF proteins [77]. It has been proposed that the origin of these genes in Drosophila is due to the recruitment for cellular functions from transposases in two independent events [77]. These PCGs can have one or more domains with high similarity and are able to encode regulatory proteins in vertebrates [78] and are associated to organ development [74, 79]. The number of BED domains is variable among ZBED genes; however, they present at least one zf-BED domain (zinc finger). The influence of zinc motifs in transcription has been largely described [80, 81], and here we present several candidates to exaptation with these motifs in genome. We found six hAT transposons (from Homo10 and Homo5 families) within genes with full BED domains, all of them within lncRNA genes. The Homo family (Homo1 to Homo11) was described in the D. mojavensis genome as having putatively active and non-active elements with length varying from ~2,500 to 3,900 bp [47]. Homo5 encodes a 479 aa transposase, has TIRs and TSD, so it is a putatively active transposon. On the other hand, Homo10 encodes a transposase of 471 aa, but it is putatively inactive because it does not have TIRs and no TSDs. The BED domains of the different genes harboring the Homo sequences share high nucleotide similarity, as identified in the ZBED genes [82] and such high similarity among lncRNA genes was observed only between lncRNA 62 genes with TE insertions harboring BED domains. This suggests an exaptation event of hAT transposons with zf-BED domains (Homo superfamily) and lncRNA genes. The high similarity may be due to functional restrictions of these genes in D. mojavensis, and also indicates that these lncRNA genes are paralogs originated from duplication events. The fact that we did not identify Homo insertions in lncRNA genes of other Drosophila species is probably due to a higher evolution rate of these genes [83]. The high similarity among all Homo insertions within lncRNA genes in D. mojavensis reinforces the functional role of TE exaptation into lncRNA genes. Our found is in agreement with recent reports that have shown that lncRNA genes have a strong association with TEs that provide a source of functional domains [37]. This association was clearly demonstrated in vertebrates, in which more than 80% of lncRNA genes harbor TE sequences [36]. Our study shows that for the D. mojavensis genome only 3.4% of lncRNA harbor TE sequences, a significantly smaller proportion than that found in vertebrates. Perhaps this is a peculiarity of these genes in Drosophila, as we also searched for TEs inserted lncRNA genes in the D. melanogaster genome and from 1,978 genes only 10.9% harbor TE insertions. Galileo and THAP domain in exons and UTRs The D. mojavensis genome harbors two Galileo insertions into exons and two in UTRs and all insertions have the THAP domain, which is characterized by a zinc- dependent DNA-binding activity [48]. The THAP domain is conserved and shared by several metazoa, and it’s related to cell-cycle regulators, pro-apoptotic factors, 63 transcriptional repressors or chromatin-associated proteins [48, 84, 85]. Two of these sequences were inserted into 5´ UTRs of PCGs and two in exons of lncRNAs. In D. melanogaster, one insertion of the DNA transposon pogo was observed in the 3’ UTR of the kinetochore Mis12-Ndc80 network component 1 (Kmn1) gene. This insertion increases the xenobiotic resistance [86]. Taking into consideration that D. mojavensis is a cactophilic species and uses columnar cacti as breeding and feeding sites [87 - 89], which presents high levels of xenobiotic components [90, 91], TE insertions at the UTR and regulatory regions could bring some advantage to D. mojavensis in this challenged environment. Although it is not possible to establish a direct relationship between Galileo insertions at the UTR of the multidrug resistance-associated protein 7 gene and the feeding habit of this species, the presence of protein domain as THAP that is DNA-binding zinc-dependent, could have an influence on transcriptional regulation. Conclusions The role of transposable elements as sources of molecular novelties is well documented and, in this study, we enlarge the examples and reinforce the importance of TEs at a functional scenario. The characterization of the mobilome of D. mojavensis at a functional perspective of transposable elements is the center of this study and shows the importance of more detailed analysis on this genome. Here we describe 14 genes candidate to harboring exapted TE sequences harboring DNA-binding sites, mostly zinc binding sites, which are inserted into exons of protein-coding genes (6) or of lncRNA genes (8). The influence of zinc motifs in transcription has been largely described [80, 81] 64 and our findings, as well as adding new evidence of TEs as source of functional domains, show that the association between specific binding sites and TEs is recurrent and maybe a main source of genomic novelties. MATERIAL AND METHODS Transposable elements annotation The TE annotation of D. mojavensis genome was assessed based upon RepeatMasker [92] annotation available at NCBI: GCF_000005175.2. All sequences that were not TEs were removed, such as simple repeats, satellite and low complexity DNA, tRNA, rRNA, unknown sequences and A, T, C, G enriched regions. To merge all LTR sequences with its respective LTR element, and also the small fractions of the same TE family that were near each other in the genome, we used the perl tool One Code to Find Them All, which parses the RepeatMasker .out file to better determine the number and positions of TE copies in the query sequence [93]. Due to RepeatMasker misclassifications caused by insufficient sequence similarity to define the correct TE family [94], we have minimized the family number super estimation by removing the species names from family classification (e.g: “Gypsy-1-Dmoj” = “Gypsy-1”). In addition, all TE insertions shorter than 100 nt were removed. 65 Identification of transposable elements associated with genes The genes annotation used to perform this analysis is from the NCBI public genome GCF_000005175.2_dmoj_caf1. We have separated each gene region into Browser Extensible Databed (BED) files. The UTRs were annotated with AGAT software (https://github.com/NBISweden/AGAT). The exon annotation is composed by the CDS regions from PCGs and exons from ncRNA genes. In order to identify TEs located at gene regions: (i) exon, (ii) 5’ and (iii) 3’ UTRs and (iv) introns, the software Bedtools [95], function intersect, was used to perform the intersection between TE insertions and gene coordinates. The parameter -f 1 was used to identify only TEs that have complete overlap with the gene regions. The TEs located at intron regions were predicted based upon the TE identification inside the gene, but not overlapping either exon or UTRs. This analysis was strand-specific; therefore, we have detected only those TE insertions that are at the same gene strand. Identification of protein domains of TEs inserted into exons All 1,593 insertions located within exons were translated into the three possible frames with EMBOSS Transeq web tool (https://www.ebi.ac.uk/Tools/st/emboss_transeq/). The translated sequences from these insertions were analyzed with Bach-CD search (https://www.ncbi.nlm.nih.gov/Structure/bwrpsb/bwrpsb.cgi) to identify their superfamily protein domains. Only the smaller score for each insertion was filtered. 66 Coding-gene harboring TE homology analysis Phylogenetic analysis To investigate the homology between one of the TEs found within exon of protein- coding genes, the DNA transposon Hopseu2, that occurred within genes presenting a ZnMc domain (Zinc-dependent metalloprotease), phylogenetic analysis was performed using the sequence of the first exon of the low choriolytic enzyme gene (LOC26528709) that have the ZnMc domain shared with the Hopseu2, and sequences such genes in other Drosophila genomes with this domain. A search for genes with this domain was made in all the Drosophila database at NCBI using BLAST and the highest hits were selected (e- value: < 3e-10; Query cover: > 50 %; identity: > 65 %). All 22 Drosophila sequences returned under the cutoff and with ZnMc domain were selected and those genomes without RepeatMasker annotation available were removed. The sequences were aligned using MAFFT [96] and a phylogenetic reconstruction was performed using nucleotide sequences from the TE sequence and those of first exon of the low choriolytic enzyme gene. The phylogeny was reconstructed using the maximum likelihood method, with 1,000 replicates and Kimura 2-parameters as substitution model in MEGAX software [97]. The accesses of gene sequences were obtained in NCBI (Additional file 1). 67 Supplementary information Additional file 1: Table S1. The NCBI access of each sequence species used to phylogenetic analyses. Additional file 2: Table S2. TE content and superfamilies percentage identified in D. mojavensis genome (193.83 Mb). TE copy number, contribution in base pairs (bp) to the genome, contributions, relative to the total TE fraction and relative to the TE copy number are given in percentages (%) of each superfamily on mobilomes. *: superfamilies not identified. Additional file 3: Figure S1. Amino acid alignment of exon 1 of the metalloprotease genes. Additional file 4: Figure S2. Phylogenetic analysis (maximum likelihood) of exon1 of genes that present ZnMc domain. The analysis was performed with 1,000 replicates in MEGAX using the Kimura 2-parameters as substitution model. It was used 22 genomes and all sequences it was obtained in NCBI. Additional file 5: Figure S3. The alignment between the Hopseu2 sequence described and identified in D. pseudoobscura genome (NC_046679.1) and Hopseu2 present in intergenic region of D. mojavensis genome. Funding This research was supported by the ANR (grant Exhyb (14-CE19-0016)) to CV and CMAC), by Brazilian agency FAPESP—Fundação de Amparo à Pesquisa do Estado de São Paulo (grant 2016/19271-2) and CNPq-Conselho Nacional de Desenvolvimento Científico e Tecnológico (grant 303455/2017-9) to CMAC, Idex Lyon and CAPES for the fellowships awarded to MCS. Availability of data and materials Data are available under request. Authors’ contributions MCS, CMAC, CV conceived the research; DSO performed bioinformatic analysis, MCS performed phylogentic analysis and analyzed the data. MCS wrote the manuscript. CMAC revised the manuscript. All authors read and approved the final manuscript. 68 Ethics approval and consent to parti