Genome-wide characterization of Mariner-like transposons and their derived MITEs in the Whitefly Bemisia tabaci (Hemiptera: Aleyrodidae).

The whitefly, Bemisia tabaci is a hemipteran pest of vegetable crops vectoring a broad category of viruses. Currently, this insect pest showed a high adaptability and resistance to almost all the chemical compounds commonly used for its control. In many cases, transposable elements (TEs) contributed to the evolution of host genomic plasticity. This study focuses on the annotation of Mariner-like elements (MLEs) and their derived Miniature Inverted repeat Transposable Elements (MITEs) in the genome of B. tabaci. Two full-length MLEs belonging to mauritiana and irritans subfamilies were detected and named Btmar1.1 and Btmar2.1, respectively. Additionally, 548 defective MLE sequences clustering mainly into 19 different Mariner lineages of mauritiana and irritans subfamilies were identified. Each subfamily showed a significant variation in MLE copy number and size. Furthermore, 71 MITEs were identified as MLEs derivatives that could be mobilized via the potentially active transposases encoded by Btmar 1.1 and Btmar2.1. The vast majority of sequences detected in the whitefly genome present unusual terminal inverted repeats (TIRs) of up to 400 bp in length. However, some exceptions are sequences without TIRs. This feature of the MLEs and their derived MITEs in B. tabaci genome that distinguishes them from all the other MLEs so far described in insects, which have TIRs size ranging from 20 to 40 bp. Overall, our study provides an overview of MLEs, especially those with large TIRs, and their related MITEs, as well as diversity of their families, which will provide a better understanding of the evolution and adaptation of the whitefly genome.

Introduction

Transposable elements (TEs) are repetitive mobile DNA sequences that can transpose from a position to another in a host genome. These elements are classified into two classes based on their transposition mechanisms (Finnegan 1989). Class I TEs, known as “copy and paste” elements, are called retrotransposons since they transpose via a reverse transcriptase that can generate cDNA from RNA intermediates. Class II TEs, also known as DNA transposons, primarily use a transposase enzyme for their excision and insertion; although some elements of this class, such as Helitron and Maverick, encode other proteins required for their transposition (Mat Razali ). A further classification for eukaryotic TEs was proposed by Wicker based on mechanistic and enzymatic differences. This classification organized the TEs in a hierarchical manner into subclasses, orders, superfamilies, and families. Since then, the TEs classification has been a challenge and several classification systems have been proposed (Curcio and Derbyshire 2003; Piégu ; Hoen ). The Tc1/Mariner superfamily is a major category of DNA transposons and is considered as one of the most widespread TE superfamilies among eukaryotes. The elements of Tc1/Mariner superfamily are classified into at least eight families according to the number of amino acid residues present between the second aspartic acid (D) and the third aspartic acid (D) or glutamic acid (E) of the transposase catalytic domain, namely: Mariner (DD34D; Jacobson ), maT (DD37D; Zhang ), VS (DD41D; Gomulski ), pogo (DDxD; Shao and Tu 2001; Dupeyron ), Tc1-like elements (TLE; DD34E; Emmons ; Sang ), TR (DD35E; Zong ), IC (DD36E; Sang ), and TRT (DD37E; Zhang ).

The members of the Mariner family (DD34D) are called Mariner-like elements (MLEs) and constitute a large family that is widespread in all organisms. MLEs are distributed into five major subfamilies mauritiana, cecropia, mellifera, elegans, and irritans, depending on their sequence similarities and phylogenetic relationships (Robertson, 2002; Bigot ). Besides, eight minor subfamilies with a limited distribution were also described, among which four have been well-described, namely vertumnana, lineata, mosellana, and drosophila, while the other four are still not well-defined (Robertson 2002; Rouault ; Filée ; Grace and Carr 2020; Dupeyron ). MLEs have a simple structure consisting of a transposase open reading frame (ORF) flanked by untranslated regions (UTRs), two terminal inverted repeats (TIRs), and TA target site duplications (TSDs; Plasterk ; Leroy ). The TIRs are generally sized from 20 to 40 bp, with the exception of the Mcmar1 element identified in the nematod Meloidogyne chitwoodi, which has large TIRs (LTIRs) sized up to 355 bp. The transposase ORF encodes a 350 amino acid polypeptide containing two major domains: (1) The N-terminal DNA-binding domain, which contains the helix-turn-helix (HTH) motif responsible for the binding of the transposase to the TIRs, and the nuclear localization sequence (NLS), which is involved in the translocation of the transposase in the nucleus (Brillet ). This domain holds also the first signature motif of the transposase, i.e., WVPREL (Augé-Gouillou ). (2) The C-terminal-catalytic domain contains the DD34D motif, meaning that the last two aspartic residues (D) are separated by 34 amino acids. The three residues are anchored to three conserved motifs, i.e., TGDE, HDNA, and YSPDLAP(x)D. This catalytic triad DD34D is required for the cleavage and the integration of the TE into the target site (Brillet ; Yuan and Wessler 2011). Only a few Mariner elements are naturally active, namely mos1 (mauritiana subfamily) detected in Drosophila mauritiana, Famar1 (mellifera subfamily) found in Forficula auricularia and Mboumar 9 (mauritiana subfamily) discovered in the ant Messor bouvieri (Jacobson ; Barry ; Munoz-Lopez and Garcia-Perez 2010; Sanllorente ). Most of the described MLEs are inactive because of the natural selection and the accumulation of mutations. Mutations can impact all parts of the elements even the ORF, making them defective and unable to encode the protein needed to process their own transposition. A fraction of these defective nonautonomous DNA transposons that kept the TIRs, TSDs, and are capable of mobilization could be considered as MITEs (Feschotte ). MITEs have the same general features of MLEs and therefore exhibit TIRs flanked by TSDs. However, they do not have an ORF coding for a functional transposase, so they rely for their transposition on the transposase of autonomous class II elements (Yang ).

TEs are major factors of genomic innovations. They confer selective advantages to the host by inducing genome sequences variation due to their movement capacity that leads to frequent DNA cleavage and rearrangement, or by being recurrent source of new functional genes and molecular regulators through a process called TEs molecular domestication (Schrader and Schmitz 2019; Jangam, ; Drezen ). For example, it was well-documented that TEs played a role in the adaptation of D. melanogaster to temperate climates (Casola ; González ). Another example in D. melanogaster was given by the insertion of a long terminal repeat of a retrotransposon into the 5ʹ end of the Cyp6g1 gene inducing its overexpression and causing the resistance of D. melanogaster to a wide range of insecticide classes (Chénais ). The cotton pest Heliothis virescens showed high levels of resistance to the Bt toxin Cry1Ac linked to the insertion of TEs into a cadherin-superfamily gene (Gahan ). In Helicoverpa armigera nine TE insertions belonging to RTE, R2, CACTA, Mariner, and hAT superfamilies hosted in exons and introns of cytochrome P450 (CyP450), glutathione S-transferase (GST), and ATP-binding cassette (ABC) transporter genes were described by Klai .

The hemipteran species Bemisia tabaci (Chen et al. 2016) is an economically important agriculture pest, which is responsible of the transmission of more than 300 plant viruses (Navas-Castillo ) and infests more than 1000 plant species including Asteraceae, Fabaceae, and Solanaceae families (Abd-Rabou and Simmons 2010). Furthermore, B. tabaci displays a high adaptability and resistance to almost all the pesticide used for its control. It resists to the organophosphates, carbamates, pyrethroids, and also to the neonicotinoids and buprofezin (Horowitz ). TEs have been described as an important source of genetic modification to acquire insecticide resistance and provide a selective advantage through their insertion sites. For example, the insertion of a MITE 0.2-kb upstream of the P450 gene CYP9M10 has been correlated to pyrethroid resistance in Culex quinquefasciatus (Itokawa ). Despite the involvement of TE in genome evolution and adaptation in many insects leading to their resistance to a wide range of insecticides, there is still no study describing TE families in the genus Bemisia.

The aim of this study was to detect Mariner-like transposons in B. tabaci genome, including their derived MITEs, and to characterize complete potently active MLEs sequences using bioinformatic tools.

Materials and methods

Supporting data

The MEAM1/B B. tabaci genome of 615 Mb available in the NCBI database (assembly ASM185493v1) was used to identify MLEs and their derived MITEs. This genome has been sequenced by Chen et al. (2016) using both Illumina short reads and PacBio long-read approaches, then assembled into 19,751 scaffolds with 3,232,964 kbp N50 length.

Identification of MLEs and their derived MITEs

The search for MLEs in B. tabaci genome was carried out according to a similarity-based method as described by Bouallègue and Xie A total of 57 nucleotide and 19 protein reference sequences from the five major Mariner subfamilies (irritans, mauritiana, mellifera, elegans, and cecropia) were used as queries (Supplementary Tables S1 and S2). The nucleotide and protein queries were used to perform BLASTN and tBLASTN searches, respectively (Altschul ). A first filtration step was conducted to exclude sequences with identity lower than 30%. The remaining sequences were extended to 1000 extra base pairs of their flanking regions in order to get full-length MLE sequences with TSDs, TIRs, UTRs, and ORF. The TIRs were manually detected by aligning the sequences with their reverse complement using multAlin algorithm (Corpet, 1988) . In order to detect the one-sided TIR sequences, we aligned them with those holding the two-side TIRs. The identified MLE sequences were next used as nucleotide queries to retrieve more sequences using BLASTN algorithm. The obtained sequences were submitted to a second filtration step to exclude redundant copies and those smaller than 500 bp (Wallau ; Berthelier ). The remaining sequences were clustered using CD-HIT with a threshold of 80% identity and a coverage of 60%. Then, their ORFs were identified using BLASTX algorithm and SIXPACK program implemented in the EMBOSS portal (https://www.ebi.ac.uk/Tools/st/emboss_sixpack). Additional searches of the HTH DNA-binding domain and NLS were performed using NPS@: Network Protein Sequence Analysis and SeqNLS (Dodd and Egan 1990; Lin and Hu 2013), respectively. The search for the catalytic domain represented by the three aspartic acid separated by 34 amino acid anchored to the TGDE, HDNA, and YSPDLA(X)D motifs was carried out manually. Finally, a rearmost filtration step was carried out to discard copies with no catalytic domain.

To retrieve MITEs derived from MLEs, the MEAM/B genome was first submitted to MITE Tracker (Crescente ) and the obtained putative MITE sequences were then used as queries to perform a BLASTN (e-value < 10−10; Lu ) against the previously identified MLEs.

MLEs clusterization

The MLE sequences were compared to each other and clustered based on a sequence identity higher than 80% and a minimum alignment coverage of 60% to their longest sequence. The obtained lineages were visualized by Cytoscape 3.7.2 software then used to perform maximum-likelihood analysis by mean of MEGAX software (Kumar ). The assignation of the MLEs lineages into their subfamilies was performed based on their comparison to the NCBI database using BLASTX algorithm, then, a maximum-likelihood phylogenetic analysis was conducted using PhyML 3.0 software package (Figure 1 ), based on the alignment of the reconstructed amino acid sequences from each lineage with known reference sequences (Supplementary Table S3), the bootstrap replicates are ranging from 0 to 1 (Guindon ). These sequences were aligned using CLUSTAL 2.1 Multiple Sequence Alignments software (Chenna ).

Figure 1

Phylogenetic relationships based on the reconstructed amino acid sequences from each lineage of the identified MLEs and the known reference sequences of mariner elements. The tree was inferred using the maximum-likelihood method with a bootstrap ranging from 0 to 1.

Transcriptome data analysis

The transcriptome data analysis was conducted in order to know whether the complete MLE sequences are expressed considering the protocol described by Pertea . In this study, we used six sequence read archive (SRA)-runs of the Whitefly B. tabaci MEAM1 transcriptome data under the bioproject accession PRJNA312467. We studied the expression of the two complete MLEs on whiteflies during the first 3 days of acquisition of Tomato yellow leaf curl virus (TYLCV).

The SRA data were first downloaded using NCBI SRA Toolkit, these SRA were mapped using HISAT2 (http://ccb.jhu.edu/software/hisat2 or http://github.com/infphilo/hisat2) to the Scaffolds where the complete MLEs were identified. The alignments are submitted to StringTie (http://ccb.jhu.edu/software/stringtie or https://github.com/gpertea/stringtie) to assemble and quantify the transcripts of each sample. Assembled transcripts were merged together by a special StringTie module, which creates a uniform set of transcripts for all samples. In order to compare the genes and transcripts with the reference annotation and reports statistics, we conducted an optional step using gffcompare (http://ccb.jhu.edu/software/stringtie/gff.shtml or http://github.com/gpertea/gffcompare). Finally, StringTie estimate transcript abundances and create table counts for Ballgown, which is an R package, used to calculate the normalized FPKM (fragments per kilobase of exon per million reads mapped) values and genes with FPKM > 1 (Macko-Podgórni ).

Results

Presence of MLEs in the genome of B. tabaci

The search for MLEs in the MEAM1/B genome led to the identification of 550 sequences. The total genomic mass of all the MLEs is 0.6 Mb, or about 0.1% of the B. tabaci genome size of 615 Mb. The NCBI database investigation using BLASTX algorithm showed that 294 of the identified elements shared the best amino acid identity, greater than 44%, with mauritiana elements while the 256 remaining elements shared the best amino acid identity, greater than 52%, with irritans elements.

To confirm these latter results, a maximum-likelihood phylogenetic analysis of the reconstructed amino acid sequences from each lineage was performed (Figure 1). Repeatedly, B. tabaci Mariner elements were found in mauritiana and irritans subfamilies indicating significant polymorphism among the different lineages of each subfamily.

Mauritiana sequences analysis

The mauritiana MLEs clusterization resulted in 6 single copies, i.e., Btmar1.1 to Btmar1.6, and 10 different lineages, i.e., Btmar1.7 to Btmar1.16 (Table 1). Clustering analysis using the Cytoscape tool showed that the majority of the mauritiana elements were clustered into the lineages Btmar1.7 (89) and Btmar1.8 (76) (Table 1). These elements are sized from 631 to 986 bp. All these copies are defective due to mutations in their ORFs, which code for a truncated transposase that have only one mutated catalytic motif. Most of the elements in the Btmar1.7 and Btmar1.8 lineages have LTIRs of 125 bp in length and flanked by TA TSDs.

Table 1

Characteristics of mariner-like transposons of the mauritiana subfamily in the whitefly genome

	Names	Number of related sequences	Sequences length (bp)	Transposase motifs	TIR length (bp)	TSD
Single copy	Btmar1.1 a	1	1,573	TGDE-HDNA-YSPDLAPAD	5ʹ:180	5ʹ: TA
					3ʹ:180	3 ʹ: TA
	Btmar1.2	1	1,325	TGDE-HDNA-YSPDLAPAD	5ʹ:217	5ʹ: TA
					3ʹ:217	3ʹ: TA
	Btmar1.3	1	1,225	YSPDLVPCD	5ʹ:287	—
					3ʹ:282
	Btmar1.4	1	1,118	HHHA-PLPLIHLTLLLV	5ʹ:308	—
					3ʹ296
	Btmar1.5	1	654	HDNA-YSPDLAPCD	5ʹ:343	—
					3ʹ:23
	Btmar1.6	1	1,188	HDNA-YSPESPDLAPCD	5ʹ:446	5ʹ: -
					3ʹ:448	3ʹ: TA
Lineages	Btmar1.7	89	631–986	FGDN	5ʹ:125	5ʹ: TA
					3ʹ:125	3ʹ: TA
	Btmar1.8	76	757–978	FGDN	5ʹ:125	5ʹ: TA
					3ʹ:125	3ʹ: TA
	Btmar1.9	34	535–657	DLAPCD	5ʹ: 130	5ʹ: TA
					3ʹ: -	3ʹ: -
	Btmar1.10	33	528–703	DLAPCD	5ʹ: -	5ʹ: -
					3ʹ:438	3ʹ: TA
	Btmar1.11	23	731–1193	HDNA-YSPDLAPCD	5ʹ:439	5ʹ: TA
					3ʹ:440	3ʹ: -
	Btmar1.12	17	642–966	FGDN	5ʹ:125	5ʹ: TA
					3ʹ:125	3ʹ: TA
	Btmar1.13	7	793–1,177	HDNA-YSPNLAPCD	5ʹ:444	5ʹ: TA
					3ʹ:227	3ʹ: -
	Btmar1.14	4	554–700	HDNA-YLSD	5ʹ:122	5ʹ: TA
					3ʹ: -	3ʹ: -
	Btmar1.15	3	1,216–1,194	HDNA-YSPDLAPCD	5ʹ:459	5ʹ: TA
					3ʹ:459	3ʹ: TA
	Btmar1.16	2	1,586–1,705	TGDE-HDNA-CSPDLAPCD	—	—

Full-length copy.

The other mauritiana copies are also defective. However, most of their transposase coding sequences retain the second and third aspartic acid of the DD34D catalytic motif. Some elements have both LTIRs retained, ranging in size from 125 to 459 bp. Whereas other elements are truncated copies with only one terminal sequence, ranging in size from 122 to 438 bp. The remainder is fragments not associated with any recognizable LTIRs, suggesting they are old and degenerated. These LTIRs could be flanked or not by TA TSDs. All the mauritiana sequences described are truncated with the notable exception for Btmar1.1, which is a full-length MLE (Table 1).

Feature analysis of Btmar1.1

The mauritiana element Btmar1.1 is placed on Scaffold 26,953, spanning from positions 130,789 to 132,365. This element revealed an ORF sized up to 1077 bp and encoding a putative protein of 358 amino acids. The comparison of this protein sequence to those in the database revealed identities with various MLE transposases, and the best identity was greater than 53% with the mauritiana element from Cryptotermes secundus (Sequence ID: PNF20350.1). As shown in Figure 2, the Btmar1.1 putative transposase displays the three catalytic motifs TGDE, HDNA, YSPDLAPAD in the C-terminal domain and a slightly different WVPREL motif (LVPRLL in Btmar1.1), a possible NLS and an HTH motif in the N-terminal domain. The Btmar1.1 transposase ORF is flanked by 5ʹ and 3ʹ UTRs, nearly perfect 180-bp LTIRs, which show 11 substitution distributed all along the LTIRs, and TA TSDs.

Figure 2

Sequence features of the full-length mauritiana element Btmar1.1 (1573 bp) and Btmar2.1 (1807 bp). The pink and blue arrows at the extremities refer to the TIRs of Btmar1.1 and Btmar2.1, respectively. TA TSDs are written with bold characters located toward the outer end of each TIRs. The transposase ORF is represented in the figure by a yellow box for Btmar1.1 and green box for Btmar2.1, respectively. The three signature motifs of the catalytic domain motifs DD34D are given within the ORF. The bipartite NLS is shown as a blue line, the HTH motif is defined by a blue diamond for the two elements.

Irritans sequences analysis

The irritans MLEs clusterization resulted in seven single copies, i.e., Btmar2.1 to Btmar2.7, and nine different lineages, i.e., Btmar2.8 to Btmar2.16 (Table 2). The majority of the irritans elements were clustered into the lineages Btmar 2.8 (120 elements) and Btmar 2.9 (87 elements) (Table 2). These elements are ranging from 605 to 990 bp in size. All these copies are defective due to mutations in their ORFs, which encode a truncated transposase that has only the two last catalytic motifs QDNA and YSPDLAPCD. Most of the elements in the Btmar2.8 and Btmar2.9 lineages have LTIRs of 210 bp in size and flanked by TA TSDs. The remaining irritans copies are also truncated and usually have either LTIRs on both sides with a size ranging from 200 to 237 bp or mutated TIR on one side with a size ranging from 27 to 197 bp. The following are exceptions, Btmar 2.7, which has no TIRs, Btmar2.3 and Btmar2.4, which have typical TIRs of 11 and 28 bp, respectively, but mutated, and especially Btmar2.1, which is a full-length MLE (Table 2).

Table 2

Characteristics of mariner-like transposons of the irritans subfamily in the whitefly genome

	Names	Number of related sequences	Sequences length (bp)	Transposase motifs	TIR length (bp)	TSD
Single copy	Btmar2.1 a	1	1,807	TGDE-HDNA-YSPDLAPSD	5ʹ:200	5ʹ:TA
					3ʹ:199	3ʹ:TA
	Btmar2.2	1	534	QDNA-YSPVLPLSD	—	—
	Btmar2.3	1	1,252	TGDE-HDNA-YSPDLAPCD	5ʹ:11	—
					3ʹ:11
	Btmar2.4	1	1,359	TGDK-HDNA-YSPDLAPSD	5ʹ:28	—
					3ʹ:28
	Btmar2.5	1	515	PDLALCD	5ʹ: -	5ʹ: -
					3ʹ:197	3ʹ:TA
	Btmar2.6	1	606	YSPMMTHRD	5ʹ: -	5ʹ: -
					3ʹ:27	3ʹ:TA
	Btmar2.7	1	605	QDNA-YSPDLAPCD	—	—
Lineages	Btmar2.8	120	605–990	QDNA-YSPDLAPCD	5ʹ:210	5ʹ:TA
					3ʹ:210	3ʹ:TA
	Btmar2.9	87	628–983	QDNA-YSPDLAPCD	5ʹ: 210	5ʹ: TA
					3ʹ: 210	3ʹ:TA
	Btmar2.10	15	702–1,017	QDNA-SYSPNLAPCD	5ʹ: 229	5ʹ: TA
					3ʹ:225	3ʹ:TA
	Btmar2.11	13	737–1,018	QDNA-YSPDLTPSD	5ʹ:237	5ʹ:TA
					3ʹ:237	3ʹ:TA
	Btmar2.12	5	743–996	QDNA-YSPDLPSD	5ʹ:233	5ʹ:TA
					3ʹ:235	3ʹ:TA
	Btmar2.13	3	986–1,000	QDNA-YSPDLPPSD	5ʹ:233	5ʹ:TA
					3ʹ:233	3ʹ:TA
	Btmar2.14	2	544–573	QDNA-YSPDLAPCD	5ʹ: -	5ʹ: -
					3ʹ: 197	3ʹ: TA
	Btmar2.15	2	527–559	YSADLAPCN	5ʹ: 197	5ʹ:TA
					3ʹ: -	3ʹ: -
	Btmar2.16	2	528–534	NSADLAPCD	5ʹ: -	5ʹ: -
					3ʹ:197	3ʹ:TA

Full-length copy.

Feature analysis of Btmar2.1

The irritans element Btmar2.1 is located on Scaffold 831, spanning from positions 3,236,801 to 3,238,607. The comparison of Btmar1.1 to Btmar2.1 did not show any significant similarity at the nucleic acid level while it revealed a highly conserved regions corresponding to the three catalytic motifs. Btmar2.1 displays an ORF up to 1092 bp in size, encoding a putative 363 amino acid protein. The comparison of this protein sequence to those in the database revealed identities with various MLE transposases and the best identity was greater than 48% with the irritans element from the coleopteran insect Anoplophora glabripennis (Sequence ID: XP_023310047.1). As shown in Figure 2, the Btmar2.1 putative transposase displays the three catalytic motifs TGDE, HDNA, and YSPDLAPSD in the C-terminal domain and a possible NLS, an HTH motif, and the first signature motif WVPREL in the N-terminal domain. The Btmar2.1 transposase ORF is flanked by 5ʹ and 3ʹ UTRs, nearly perfect large 200-bp TIRs, and TA TSDs. The Btmar2.1 LTIRs exhibit fewer substitutions than those of Btmar1.1, i.e., only four at positions 18, 179, 193, and 198, and a single deletion of one nucleotide at the 36th position of the 3ʹ LTIR.

MITEs deriving from MLE: detection and sequence analysis

The search for MITEs in the MEAM1/B genome led to the identification of 8024 MITE families. Among these families, 1030 MITE families belong to the Tc1/mariner superfamily from which we identified 71 MITE families related to MLE sequences. These 71 MITE sequences showed 12,605 copies, TA TSD patterns, and TIRs. BLASTN (e-value < 10−10) sequences analysis revealed that a total of 35 MITEs families derived from mauritiana elements and 36 MITEs derived from irritans elements. The MITE sequences showed different groups based on their TIR sequences (Figure 3). Most of the MITEs derived from mauritiana MLEs originated from mutated MLEs sequences belonging to Btmar1.7, Btmar1.8, Btmar1.12, and Btmar1.15 lineages. These elements are ranging from 128 to 790 bp in size with TIRs ranging from 4 to 314 bp in size. However, most of the MITEs derived from irritans MLEs originated from the full-length element Btmar2.1. These elements are ranging from 179 to 796 bp in size with TIRs of 8 to 282 bp.

Figure 3

Characteristics of MLE-derived MITEs in the whitefly genome. (A) MITEs derived from the mauritiana subfamily. (B) MITEs derived from the irritans subfamily. Light blue rectangles refer to mauritiana MLEs while light blue diamonds refer to irritans MLEs. MITEs are represented by ellipses. Green ellipses mean that MITEs have four possible counterparts but they are derived from one of them. Pink ellipses mean that MITEs have three possible MLEs origins, orange and yellow ellipses indicate that MITEs have two or one possible MLEs origins.

Sequences alignments revealed that in all the detected MITEs the first 10 base pairs in TIRs are conserved with the exception of 18 copies that have only 4–8 bp conserved in TIRs, while 20 elements have conserved LTIRs with a size exceeding 100 bp and up to 314 bp. All the described MITEs exhibited mutations such as deletions and substitutions targeting the middle of the sequences and extended to the extremities in highly mutated MITEs, which led to the generation of the MITEs with imperfectly conserved TIRs of about 4–20 bp.

Transcriptome analysis of Btmar1.1 and Btmar2.1

Btmar1.1 is inserted in the 9th intronic region of the cytosolic carboxypeptidase2-like gene (LOC109039198), which codes for the CCP2 enzyme involved in the deglutamylation reaction on tubulin, leading to Δ2-tubulin (Tort ). We tested the expression of this gene in order to see whether Btmar1.1 is expressed.

The TYLCV transcriptome analysis showed that the LOC109039198 gene has an FPKM higher than one for the dix SRA-runs (Supplementary Table S4), which means that this gene is expressed. Accordingly, the LOC109039198 gene, and then Btmar1.1, is also expressed.

The FPKM for Btmar2.1 is also higher than one in the same transcriptome conditions (Supplementary Table S4), which means that Btmar2.1 is expressed.

Discussion

In this study, we have shown that the MEAM/B genome exhibit a total of 550 MLEs belonging to both mauritiana (Btmar1.1 to Btmar1.16) and irritans (Btmar2.1 to Btmar2.16) Mariner subfamilies. The majority of Mariner sequences are mainly represented by four lineages (Btmar1.7, Btmar1.8, Btmar2.8, and Btmar2.9) comprising 67% of the identified Mariner elements. Although high copy number lineages have been described in some species such as the Tricladida Girardia tigrina, the Hemiptera Rhodnius prolixus (Garcia-Fernandez ; Filée ; Fernández-Medina ), and the Drosophila genus genomes (Wallau ), MLEs are usually described as a low copy number family of TEs and most of the other genomes display less than 50 copies per genome, hence, the B. tabaci genome seems to be rather permissive for MLEs expansion.

As previously described, the typical size of MLEs is 1300 bp (Robertson and Martos 1997), which is consistent with the size of most of the identified MLEs described in this study. Among the 550 identified elements, only 4 MLEs sequences have an unusual length, which exceeded the known size by 277–507 bp. However, analysis of the TIRs has shown that their length is extended (up to 459 bp) in most sequences. This particularity is usually observed among the closely related TLE transposons (Robertson, 1995; Shao and Tu 2001; Ruvolo ). Although TLE have variable TIRs size, ranging from 24 to 756 bp, most of them are over 50 bp in length (Ruvolo ; Jehle ). In addition, MLEs with LTIRs have been described previously in the genome of the phytoparasitic nematode M. chitwoodi, and these elements, namely Mcmar1-1 and Mcmar1-2, present perfectly conserved LTIRs of 355 bp (Leroy ).

In this study, most of the MLEs sequences are degenerated copies with in-frame stop codons and frameshift, or indel mutations. As reported previously, these mutations are sufficient to abolish the catalytic activity or prevent the synthesis of the putative transposase and make the MLEs inactive or nonautonomous (Lohe ). Importantly, our results have also highlighted two complete copies (from TSD to TSD) with putative conserved transposases, which make them probably capable of their own transposition and the transposition of close or unrelated nonautonomous elements via trans-complementation (Lohe ).

As previously reported, MITEs can be derived from autonomous or nonautonomous elements due to mechanisms such as abortive gap repair (de Ortiz ; Fattash ). However, the lost or the inactivation of their original autonomous elements does not make them imperatively inactive, a few MITE sequences can be cross-mobilized using nearly noncognate autonomous elements machinery (Fattash ). For example, the loss of the original autonomous element Ping did not stop the activity of the MITE element mPing and this element can be cross-mobilized via the transposase of other elements such as Pong (Yang ). Moreover, a mariner derived Stowaway-like MITE element, named MiHsmar1 has been mobilized using the reconstructed Hsmar1 transposase (Miskey ). For most of the identified MITEs-MLE sequences, we have been able to retrieve relatively intact TIR sequences and TA TSDs needed for the transposition. The TIRs presented high similarities with their related mariner subfamily. These findings support that in B. tabaci genome, MITEs can be mobilized through the transposase of mariner lineages that share almost identical TIR patterns.

MITEs have been widely described in plant genomes as small elements shorter than 800 bp and with a high copy number reaching thousands of copies (Guermonprez ; Crescente ). MITEs have also been characterized in insects. For example, the genome of Drosophila sechellia harbors only 46 MITE sequences (Dias and Carareto 2011). However, in the genome of Aedes aegypti the MITE copy numbers range from 2100 to 3000 per haploid genome (Tu 1997). In a more recent study, Han identified 6012 MITE families from the genome of 98 insect species and constructed the first insect MITE database (iMITEdb). In this database, the number of the MITE families belonging to the Tc1/mariner superfamily is 1698 families from 90 insect species genomes, among which 192 families from 8 Hemiptera genomes. Based on the data obtained in this study, the relative abundance of MITEs in the B. tabaci genome is high compared to other insect genomes since the number of the identified MITE related to MLE sequences exceeds the total MITE copies number described in the majority of the insects. Although, most of the identified MITEs did not originate from autonomous elements, MITEs related to full-length TEs have been described in some species such as D. sechellia (Dias and Carareto 2011). Here, we also highlight four copies of MITEs related to the full-length element mauritiana and 20 related to the irritans full-length element.

Conclusion

This work represents the first in silico characterization and classification of MLEs and their related MITEs in the genome of B. tabaci. The results highlight the presence of LTIR MLEs and their related MITEs, as well as the presence of two full-length MLEs that may be involved in the transposition of nonautonomous elements. These results provide new data for further analysis to understand the involvement of TEs in the adaptive evolution of the whitefly.

Data availability

Supplementary materials are available at GSA figshare portal: https://doi.org/10.25387/g3.15001164. Supplementary Tables S1 and S2 contain nucleic acid and transposase MLEs queries for BLASTN and tBLASTN searches, respectively. Supplementary Table S3 contains phylogenetic tree references characteristics. Supplementary Table S4 contains SRA data analysis of the expression of the two complete MLEs in Bemisia tabaci. All identified MLE and MITE sequences are available in text files S1 and S2, respectively. The authors affirm that all data necessary for confirming the conclusions of the article are present within the article, figures, and tables.

Funding

This work was financially supported by the « PHC Utique » program of the French Ministry of Foreign Affairs and Ministry of Higher Education, Research and Innovation and the Tunisian Ministry of higher education and scientific research in the CMCU project number « 19G0901 ».

Conflicts of interest

The authors declare that there is no conflict of interest.