A miniature inverted-repeat transposable element, AddIn-MITE, located inside a WD40 gene is conserved in Andropogoneae grasses

MITE identification

The identification of MITE families in 674 Mb of sugarcane genomic assembly (Grativol et al., 2014) was performed through a computational approach for de novo identification of these elements - MITE-Hunter pipeline (Han & Wessler, 2010). The pipeline of the MITE-Hunter has four main steps: (i) identification of MITE candidates; (ii) filter out of false positives through a pairwise sequence alignment; (iii) selection of MITE examples; (iv) filter out of false positives through a multiple sequence alignment; (v) grouping of predicted MITE consensus sequences into families.

The MITE-Hunter was run with default parameters with the sugarcane genome sequence file (sugarcaneMFscaffolds-1.0.fa.gz) downloaded from http://lbmp.bioqmed.ufrj.br/genome/download_sequence. The predicted sugarcane MITE sequences were manually inspected for the presence of TSD (target site duplications). Only MITEs with identified TSDs were selected for further analysis. The fasta format sequences of identified MITEs are available at Data S1.

Identification and annotation of genes associated with MITEs in sugarcane

In total, 145 sugarcane scaffolds containing MITE sequences with TSDs were aligned against Sorghum bicolor version 3.1 coding sequences (CDS) by BLASTN (Altschul et al., 1990) with the E-value cutoff of 10⁻¹⁰. Only the best hits were considered for identification of genes nearby MITEs. The sorghum CDSs and their annotation were downloaded from JGI Phytozome 12.

To confirm the annotation of exons regions in the scaffold5050—size5174 containing the gene with a WD40 domain, we used the raw sugarcane RNA-seq reads of four control samples (Vargas et al., 2014; Santa Brigida et al., 2016) downloaded from NCBI SRA under accessions (SRX374577, SRX374581, SRX603441 and SRX603445). The quality of the libraries was evaluated using FASTX Toolkit (http://hannonlab.cshl.edu/fastx_toolkit/) and reads having base quality greater or equal to 20 (Q20) was used for additional analyzes. Filtered RNA-seq reads were mapped using Bowtie2/2.1.0 (Langmead & Salzberg, 2012) onto the WD40 gene with default settings.

WD40 annotation and phylogenetic tree

For the construction of the phylogenetic tree, WD40 protein sequences from sorghum and maize were obtained in Phytozome 12. We used the annotation file of sorghum and maize genes to select only those genes annotated with WD-40 or WD repeat domain. This WD-40 domain (PF00400) was confirmed in all protein sequences from selected genes using Pfam (https://pfam.xfam.org/). Additionally, we performed a BLASTN search on a sugarcane transcriptome (TR7) (Santa Brigida et al., 2016) using the sugarcane Scaffold5050 annotated as WD40 protein. One transcript from the Locus_26251 was selected. The protein sequence predicted from this sugarcane transcript were obtained from Expasy translate tool. All selected protein sequences were aligned using MUSCLE (Edgar, 2004).

The evolutionary history was inferred by using the Maximum Likelihood method based on the Whelan and Goldman model (Wheland & Goldman, 2001). Initial tree(s) for the heuristic search were obtained automatically by applying Neighbor-Join and BioNJ algorithms to a matrix of pairwise distances estimated using estimated using a JTT model approach, and then selecting the topology with superior log-likelihood value. The rate variation model allowed for some sites to be evolutionarily invariable ([+I ], 3,37% sites). The analysis was performed by 1,000 generations and a consensus tree was generated to assign an a posteriori probability to each node using the 1,000 trees sampled. The tree was drawn to scale, with branch lengths measured in the number of substitutions per site. The analysis involved 12 amino acid sequences. All positions containing gaps and missing data were eliminated. There were a total of 274 positions in the final dataset. Evolutionary analyses were conducted in MEGA7 (Kumar, Stecher & Tamura, 2016).

The gene structure and synteny analysis of Sobic.001G093800 and GRMZM2G056645 were performed at Gene Structure Display Server v2.0 and Plant Genome Duplication Database (http://chibba.agtec.uga.edu/duplication/), respectively.

Validation of AddIn-MITE position at sugarcane genome

Following the Kuijper’s leaf numbering system for sugarcane plants (Cheavegatti-Gianotto et al., 2011), leaf -2 tissue of each sugarcane wild species and cultivars were collected for DNA extraction. Young leaves were collected from plants of S. Spontaneum clone SES205A, S. officinarum clone 82–72, S. robustum clone Molokai5009, S. barberi clone Khagziand, S. sinense clone Chukche maintained in the germplasm collection of Instituto Agronômico de Campinas, Ribeirão Preto and Monsanto Breeding Station, Maceió. Samples from SP70-1143 and other cultivars were collected from our labs’ collection. Total genomic DNA was extracted from leaves using the CTAB method (Doyle & Doyle, 1987) with minor modifications. The quality and concentration of DNA were estimated using Thermo Scientific NanoDrop™ 2000c Spectrophotometer and then the integrity was verified by electrophoresis on a 1% agarose gel.

PCR amplification was performed using 5ng of DNA from each sample. Specific primers for the validation of MITEs location were designed using the Primer3 program and are listed in Table S1. The amplification reactions were performed with a final volume of 25 µL. In each reaction, 50 mM KCl, 10 mM Tris-HCl (pH 8.0), 1 mM MgCl2, 0.2 mm dNTPs (Promega^®, Madison, WI, USA), 1U Taq DNA polymerase (New England Biolabs) and 200 nM of the primer were used. The PCR program started with an initial denaturation at 94 °C for 5 min followed by 12 touchdown cycles of 30 s at 94 °C; 45 s of annealing temperature, starting at 52 °C and decreasing 0.5 °C per cycle; and 2 min at 72 °C. Another 30 cycles with annealing temperature fixed at 46 °C and 10 min at 72 °C were used as a final step. Amplified products were observed by electrophoresis in 0.5 × TAE at 100 V using 2% agarose gel with a 100-bp ladder as standard (Promega).

Conservation analysis of sugarcane AddIn-MITE

Estimation of AddIn-MITE copy number and conservation among monocots’ genomes were performed through TARGeT: Tree Analysis of Related Genes and Transposons (Han, Burnette & Wessler, 2009). The BLASTN tool was run with sugarcane MITE sequences against Brachypodium dystachyon, Oryza sativa, Panicum virgatum, Setaria italica, Sorghum bicolor and Zea mays genomes. The coverage of the sugarcane AddIn-MITE and number of hits on each genome were plotted. To test the co-localization of the AddIn-MITE with genes in the various monocots genomes, the AddIn-MITE sequence was used as the query on NCBI/BLASTN against non-redundant (nr) sequences. All genes containing an AddIn-MITE inserted in genic regions or flanking sequences were selected and the location of the AddIn-MITE was analyzed.

Analysis of small-RNAs derived from AddIn-MITE

Small RNA libraries from leaves of sugarcane (GSM1040783), sorghum (GSM803128) and maize (GSM433620) were downloaded from NCBI’s Gene Expression Omnibus. The sRNA reads from each species were aligned against the MF scaffolds from sugarcane and the WD40 genic region associated with AddIn-MITE from maize and sorghum, using the sequence alignment tool from UEA sRNA toolkit-Plant Standalone version. The alignments were performed with non-redundant sRNAs, zero mismatches allowed and visualized at the same toolkit. Next, the aligned sRNAs were distributed by size. The hairpin structures of AddIn-MITE on the three species were obtained on RNAFold Web Server.

To measure the abundance of AddIn-MITE-derived sRNAs generated in sugarcane plants subject to pathogen infection and salt stress, the sRNA libraries under accessions codes GSM1350704, GSM1350705, GSM1040783, and GSM1040786 were downloaded from NCBI’s Gene Expression Omnibus. The alignments of WD40 genic region and sRNAs were performed as described above. The venn diagram was constructed with the sRNAs sequences mapped on AddIn-MITE region.

Identification of MITEs in sugarcane genomic assembly

The search for MITEs in the sugarcane genomic assembly of methyl-filtered reads (MF scaffolds) was performed following the MITE-Hunter pipeline, which was run using 1,109,444 sugarcane MF scaffolds (Grativol et al., 2014). Five major steps constitute the pipeline and were used to obtain 20 TE consensus sequences (Fig. 1A). An additional pipeline was applied to select TEs with TIR-like structures flanked by TSDs, which have 2–10 bp. After that, eight valid TE exemplars were selected from the 20 TE consensus sequences, which distributed among 145 sugarcane MF scaffolds (Table 1). The most common TSDs size was 3 bp for each candidate. Figure 1B shows the presence of 2 bp TSDs in the boundaries of one of the identified MITEs. The size of identified ranged from 106 to 155 bp. The sequences of sugarcane MITEs were compared with the structures of plant MITEs already described in the literature. We found that sugarcane_1_16642 has a similar structure of Stowaway MITE, as described by Bureau & Wessler (1994). Six sugarcane MITEs were similar to Tourist superfamily (Bureau & Wessler, 1992; Zhang et al., 2004).

Gene-associated MITEs

Evidence that MITEs could be associated with genic regions in the sugarcane genome were indicated by the low methylation level of this class of transposon when compared to retrotransposons (Grativol et al., 2014). In the Arabidopsis genome, MITEs are distributed throughout euchromatic regions (Wright, Agrawal & Bureau, 2003). Based on this, we have searched for genes located in neighboring regions of sugarcane MITEs. The 145 scaffolds containing eight distinct MITE families with TSDs were submitted to BLASTN against sorghum CDS to predict exonic regions on sugarcane genomic sequences. Sixty-six scaffolds showed significant alignments with sorghum CDSs (Fig. 1A). Detailed analysis of MITE insertions through the start-end position of sorghum CDSs showed that the majority of them are localized at intronic regions of genes. Only three MITEs are localized upstream and five are downstream of the CDSs (Table 2). MITEs were located in genes involved with plant development (e.g., WD40-REPEAT FAMILY PROTEIN, REPLICATION FACTOR C1 and VACUOLELESS1); hormone response (e.g., AUXIN-RESPONSIVE PROTEIN IAA1); cell wall formation (e.g., CALLOSE SYNTHASE 8 and FASCICLIN-LIKE ARABINOGALACTAN PROTEIN 17); stress response (e.g., CALCIUM-DEPENDENT PROTEIN KINASE 4, LIPID-TRANSFER PROTEIN 1, TOBAMOVIRUS MULTIPLICATION PROTEIN 3 and VASCULAR PLANT ONE ZINC FINGER PROTEIN); and epigenetic modifications (e.g., HAC13 HISTONE ACETYLTRANSFERASE and O-METHYLTRANSFERASE FAMILY 2 PROTEIN) (Table 2). The MITE family with greater copy number within genic regions of sugarcane was the sugarcane_1_16642, which we call AddIn-MITE (a piece of sequence that can be added to a region to give extra features or functions) (Table 2).

Once we found that the AddIn-MITE was associated with sugarcane genic regions, we extended the analysis to other monocot genomes. First, the conservation of the AddIn-MITE sequence in the other six-monocot genomes was analyzed through TARGeT. The analysis of the number of hits and the coverage showed that the AddIn-MITE is conserved in closely related species of Panicoideae sub-family, such as Panicum virgatum, Setaria italica, Sorghum bicolor and Zea mays (Fig. 2). In addition, BLASTN of the AddIn-MITE sequence against the non-redundant database showed that this MITE is associated with genic regions also in these four genomes (Table S2). In Z. mays, S. bicolor, S. italica and S. viridis, AddIn-MITE have similar localization inside genes. The analysis of AddIn-MITE position on those genes revealed that it is located in intronic regions (Table S2). However, no conservation was observed in species outside the Panicoideae sub-family such as Oryza sativa and Brachypodium distachyon (Figs. 2A, 2B)

AddIn-MITE conserved position in WD40 genes of sugarcane, sorghum, and maize

Because sugarcane commercial cultivars have complex, aneuploidy genome, with a history of several rounds of crossings (Grivet & Arruda, 2002), we investigated whether the AddIn-MITE location nearby WD40 gene was conserved in Saccharum wild species and hybrids. Figure 3A shows the location of the AddIn-MITE close to the WD40 exon 14 at the scaffold5050 and the position where the primers were designed for PCR. Three primers were designed: (i) one inside the AddIn-MITE; (ii) one inside the CDS and (iii) one outside both the AddIn-MITE and CDS regions (Fig. 3A). PCR reactions with gDNA from the SP70-1143 cultivar validated the distance among the AddIn-MITE and CDS region in the scaffold (Fig. 3B).

The position of AddIn-MITE nearby the WD40 exon was also verified in five wild sugarcane species and seven hybrid cultivars (Fig. 3C). All wild species and cultivars showed the amplified products with expected sizes for the distance between AddIn-MITE and WD40 exon (Fig. S1). Additionally, we confirmed the position of the AddIn-MITE inside the CALS8 gene in the wild species and cultivars (Fig. S2).

To confirm the close relationship between AddIn-MITE and genes, we analyzed the WD40 gene annotated on the sugarcane scaffold5050 (Table 2). First, we selected all proteins from the largest transcripts annotated as WD40 repeat family protein from sorghum and maize on Phytozome 12. The presence of the WD40 motif (PF00400) was verified in each protein sequence using Pfam. This resulted in 169 WD40 proteins in sorghum and 235 in maize. From these, 113 and 176 respectively showed only the WD40 motif. Next, we performed a BLASTN search on the sugarcane assembled transcriptome (TR7) (Santa Brigida et al., 2016) using the scaffold5050 to obtain the possible transcript from this region. We found one locus (Locus_26521) with zero e-value. The Pfam analysis of the predicted protein from the locus_26521 showed the presence of three domains: WD-40 (PF0400), BEACH (PF02138) and PH_BEACH (PF14844). These domains were characterized in the subfamily C of WD40 proteins in Setaria italica (Mishra et al., 2014). Considering the subfamily C of WD40 in S. italica, we selected all sorghum and maize proteins containing WD-40 and BEACH domains for a phylogenetic analysis. A phylogenetic tree was constructed with these WD40 proteins using the maximum-likelihood method (Fig. 4). Three different subgroups were found on the WD40 subfamily C. The proteins (Sobic.001G093800.1, GRMZM2G056645_T01, and Locus_26521_transcript_1) formed one subgroup, which confirmed the annotation of the scaffold5050 (Table 2).

The gene structure of Sobic.001G093800 and GRMZM2G056645 was compared with GSDS web tool (Fig. 5). The WD40 gene structure was similar between sorghum and maize. We also compared the Locus_26521_transcript_1 and Scaffold5050 predicted exons with those sorghum and maize genes (Fig. 5). A second annotation verification in WD40 exons regions at the scaffold5050 using RNA-seq reads was performed. The transcriptome libraries from sugarcane plants grown in control hydroponics and in vitro conditions were aligned against the scaffold5050 containing the WD40 exons and AddIn-MITE. The Locus_26251 assembled from RNA-seq data of sugarcane (Santa Brigida et al., 2016) and the peaks of transcription at the two exons regions (Fig. 4), confirmed the position of the WD40 exons on the sugarcane scaffold5050.

We also verified that AddIn-MITE was located close to WD40 exons in both sorghum and maize (Fig. 5). Then, we analyzed the flanking sequences and the orientation of AddIn-MITE. The orientation of AddIn-MITE was the same in the three species and their flanking sequences showed great conservation (Fig. 6A). The analysis of synteny on Plant Genome Duplication Database (PGDD) reveals a huge block on chromosome 1 of maize and sorghum, which include the WD40 genes of both species (Fig. 6B). These results suggest that the AddIn-MITE inserted close to WD40 exons before the divergence of sugarcane, sorghum, and maize.

Profile of AddIn-MITE-related sRNAs on sugarcane WD40 gene

Since the first report by Llave and coworkers in Llave et al. (2002), the knowledge about the roles of sRNA regulating gene expression in plants has been expanded enormously. In recent years, the association between repetitive genomic regions and sRNA has been described. TE-derived sRNA has an important role in the feedback silencing of the active TEs (Ito, 2013). To evaluate if the AddIn-MITE located inside WD40 gene could generate MITE-derived sRNAs, we used sRNAs libraries constructed from sugarcane leaves and mapped them to the sugarcane genomic assembly (MF scaffolds) with zero mismatches. The sRNA mapping showed that on the scaffold5050 the majority of the mapped sequences are derived from the AddIn-MITE region (Fig. 7). The distribution of sRNA size showed that the sRNAs with 21-nt were the most abundant among the mapped sequences (Fig. 7). The WD40 genes containing the AddIn-MITE in sorghum and maize were also checked for the alignment of sRNAs from leaf libraries. The mapping on sorghum showed very few sRNAs on AddIn-MITE region, most of them with 21-nt in size (Fig. S3A). Curiously, the WD40 genes from maize did not show sRNAs mapped in the AddIn-MITE region with zero mismatches (Fig. S3B).

Since the majority of canonical plant miRNA has 21-nt in length and are produced from a hairpin precursor (Lee et al., 2015), the 21-nt sequences mapping on the sugarcane scaffold5050 could be the initial step for the emergence of a new miRNA precursor. To evaluate this, we compared the hairpin structure from each AddIn-MITE region of maize, sorghum, and sugarcane (Fig. 8A). The hairpin structure in the AddIn-MITE region was found on sugarcane WD40 gene and showed the lowest Minimum Free Energy (MFE = −62.9). The sugarcane AddIn-MITE on WD40 gene also showed a pattern of 21-nt small RNA mapping on its 3′ end (Fig. 8B), which would suggest that this region is a “proto-miRNA” locus in sugarcane. This seems to be corroborated by the profile of sRNAs mapped in the AddIn-MITE region on WD40 gene. Seventeen sRNAs from sugarcane plants subjected to pathogen infection and to salt stress mapped to AddIn-MITE (Fig. S4A). From these, nine have 21-nt in length (Fig. S4).