First Report of Aleurocanthus spiniferus on Ailanthus altissima: Profiling of the Insect Microbiome and MicroRNAs

We report the first occurrence of the orange spiny whitefly (Aleurocanthus spiniferus; OSW) on the tree of heaven (Ailanthus altissima) in Bari, Apulia region, Italy. After our first observation in 2016, the infestation recurred regularly during the following years and expanded to the neighboring trees. Since then, we have also found the insect on numerous patches of the tree of heaven and other plant species in the Bari province. Nevertheless, the tree of heaven was not particularly threatened by the insect, so that a possible contribution by OSW for the control of such an invasive plant cannot be hypothesized hitherto. This work was also aimed at profiling the microbiome of OSW feeding on A. altissima. For this purpose, we used the denaturing gradient gel electrophoresis (DGGE) and the deep sequencing of small RNAs (sRNAs). Both techniques unveiled the presence of “Candidatus Portiera” (primary endosymbiont), Wolbachia sp. and Rickettsia sp., endosymbionts already reported for other Aleyrodidae. Deep sequencing data were analyzed by four computational pipelines in order to understand the reliability of the detection of fungi, bacteria, and viruses: Kraken, Kaiju, Velvet, and VelvetOptimiser. Some contigs assembled by Velvet or VelvetOptimiser were associated with insects, but not necessarily in the Aleurocanthus genus or Aleyrodidae family, suggesting the non-specificity of sRNAs or possible traces of parasitoids in the sample (e.g., Eretmocerus sp.). Finally, deep sequencing data were used to describe the microtranscriptome of OSW: 56 canonical and at least four high-confidence novel microRNAs (miRNAs) were identified. The overall miRNA abundance in OSW was in agreement with previous works on Bemisia tabaci, and bantam-3p, miR-276a-3p, miR-317-3p, miR-750-3p, and mir-8-3p were the most represented miRNAs.


Introduction
The orange spiny whitefly (OSW), Aleurocanthus spiniferus (Quaintance, 1903) originated in tropical Asia and has spread to Africa, Australia, and the Pacific Islands. In Europe, it was reported for the first time in Italy [1], and it is spreading all over the Mediterranean Countries [2,3]. Aleurocanthus Quaintance and Baker is a paleotropical genus, currently including about 80 described species [4].
In this work, we documented the first occurrence of OSW on A. altissima and used sRNA deep sequencing to obtain the microbiome profile and sRNA transcriptome.

Insect Identification
OSW puparia were collected yearly since July 2016 from diverse host plants, including the tree of heaven, in the Bari province. More than one hundred puparia were slide-mounted after discarding the layering previous instar exuvia [29]. The slide mounting procedure (https://zenodo.org/record/3471649#.Xk1cr2hKhPY) made the puparia cuticle clear enough to discriminate relevant details. Thus, the mounted puparia were scrutinized and identified using a compound light microscope to study taxonomic characters [1,4].

Denaturing Gradient Gel Electrophoresis (DGGE)
Before DNA extraction, 10 single specimens of OSW infesting A. altissima were surface sterilized in 70% ethanol and rinsed three times in sterile distilled water. DNA was extracted with a slight modification of the protocol described by Gebiola et al. [30], where specimens were homogenized with a sterile pestle directly in the extraction buffer (Chelex 100 + pK). The V3 hypervariable region of 16S rDNA was amplified by PCR [31,32] and the products were analyzed by a slight modification of the DGGE technique described by Muyzer et al. [33] using a 40%-65% denaturing gradient (90 V for 20 h). Known endosymbionts such as Wolbachia Hertig from Encarsia formosa Gahan

Insect Identification
OSW puparia were collected yearly since July 2016 from diverse host plants, including the tree of heaven, in the Bari province. More than one hundred puparia were slide-mounted after discarding the layering previous instar exuvia [29]. The slide mounting procedure (https://zenodo.org/record/ 3471649#.Xk1cr2hKhPY) made the puparia cuticle clear enough to discriminate relevant details. Thus, the mounted puparia were scrutinized and identified using a compound light microscope to study taxonomic characters [1,4].

Denaturing Gradient Gel Electrophoresis (DGGE)
Before DNA extraction, 10 single specimens of OSW infesting A. altissima were surface sterilized in 70% ethanol and rinsed three times in sterile distilled water. DNA was extracted with a slight modification of the protocol described by Gebiola et al. [30], where specimens were homogenized with a sterile pestle directly in the extraction buffer (Chelex 100 + pK). The V3 hypervariable region of 16S rDNA was amplified by PCR [31,32] and the products were analyzed by a slight modification of the DGGE technique described by Muyzer et al. [33] using a 40%-65% denaturing gradient (90 V for 20 h). were excised, eluted, used as templates for PCR with the primers 341f/518r, and directly sequenced. Obtained sequences were blasted against the GenBank database.

Deep Sequencing of Small RNAs
Total RNA was extracted from 100 mg of OSW puparia using TRIzol™ (Thermo Fisher Scientific, Bedford, MA, USA), according to the manufacturer's instructions, subjected to DNase digestion using TURBO DNA-free™ Kit (Thermo Fisher Scientific, Bedford, MA, USA), and suspended in nuclease-free water. OSW puparia were rinsed three times in sterile distilled water before RNA extraction. One microgram RNA was sent to IGA Technology Services (Udine, Italy) for library preparation and sequencing. The TruSeq Small RNA Library Prep Kit (Illumina, San Diego, CA, USA) was used for library construction, and 50 bp reads were sequenced on HiSeq 2500 system (Illumina, San Diego, CA, USA) in single-end mode.
The sequencing raw reads were subjected to adapter removal, quality control and filtering using FASTX-toolkit version 0.1 [34], FASTQC version 0.72+ (http://www.bioinformatics.babraham.ac.uk/ projects/fastqc/), and Filter FASTQ tool version 1.1.1 [35], respectively, in order to obtain high-quality reads of 15-35 nt to be used for the downstream analysis. These steps were performed within a locally-installed Galaxy version 19.01 [36].

Microbiomics
High-quality sequencing reads of 15-35 nt were used to infer the complete microbiome of OSW, i.e., the presence of reads associated with insects, fungi, bacteria, and viruses. Since the taxonomic assignment of sequencing reads might be affected by the bioinformatics approach, including the reference database, we used four computational pipelines and compared their results: Kraken, Kaiju, Velvet, and VelvetOptimiser (the latter two tools were followed by a blastn search). Kraken is a taxonomic sequence classifier that assigns taxonomic labels to short DNA reads [37]. We subjected the high-quality sequencing reads to version 1.2.4 of this tool within the public Galaxy server (https://usegalaxy.org/), where it utilizes databases for bacteria, fungi or plasmids. Kaiju is a web-based program (http://kaiju.binf.ku.dk/) for taxonomic classification of sequencing reads from metagenomic whole genome sequencing or metatranscriptomics experiments [38]. We used version 1.7.2. with default parameters and the NCBI nr database including non-redundant proteins of bacteria, archaea, viruses, fungi, and microbial eukaryotes. Velvet 1.2.10.0 was used to assemble reads into contigs with a hash length of 21 [39]. The hash length value affects the number and length of contigs finally obtained: smaller hash lengths yield more contigs, but shorter in size, whereas higher hash lengths yield fewer contigs, but longer. VelvetOptimiser assembles reads using a range of hash length values in order to maximize the contig length, contig number, N50 value or other parameters as specified by the user [40]. We used this tool with a range of hash length between 19 and 31. Contigs assembled with Velvet or VelvetOptimiser were then used as queries for the search against NCBI BLAST nucleotide collection database (nr/nt) using the blastn algorithm on the web service. Only BLAST hits with more than 95% sequence similarity and e-value lower than 10 −5 were retained.
The outputs of these four pipelines were processed with MEGAN6 for the graphical representations [41].

Classification of Small RNAs
In order to have an overview of the sRNA features in OSW, we carried out several analyses. Length distribution of reads was extracted from the FASTQC output, imported into Microsoft Excel Professional 2016 and used to generate histograms. This procedure was done both for redundant and non-redundant reads, the latter obtained using Collapse sequences tool version 1.0.1 [34] within the locally-installed Galaxy. Sequencing reads were used as queries to search sequence similarities in the Rfam database, separately by RNA categories such as rRNAs, miRNAs, several classes of repetitive elements, transcripts and other Rfam families. Alignment to the B. tabaci genome using bowtie 2 with default parameters [42] was used to determine the number of reads matching the genome. B. tabaci isolate MEAM1 (assembly ASM185493v1; 19,751 scaffolds) was used because it is a species in the same family (Aleyrodidae) of OSW, which genome is not yet available. The remaining reads were considered as unassigned RNA. The number of reads associated with these features was used to generate pie charts by Microsoft Excel Professional 2016.

Identification of Known and Novel MicroRNAs Including IsomiRs
miRDeep2 was used to identify miRNAs in OSW [43]. As inputs to this tool, we submitted the sequencing reads and the genome of B. tabaci isolate MEAM1 (assembly ASM185493v1), without known miRNA and precursor sequences. MiRDeep2 predicts precursors and miRNAs in a genome and associates them with a 0-10 score and a randfold significance. No arbitrary threshold was used to select sequences, and all predicted miRNAs were retained. MicroRNA precursors available in miRBase (release 22.1, October, 2018) [44] were mapped to the B. tabaci genome and compared to the miRDeep2 output in order to discriminate known and novel microRNA loci. Also, mature miRNAs available in the same database were aligned with the sequences predicted by miRDeep2 using the blastn algorithm. We considered as canonical miRNAs of OSW the reads with a perfect match or sharing up to two mismatches with the sequences in the database, regardless of their length. According to the miRBase guidelines on miRNA annotation [45,46] and with the concept of miRNA isoforms (isomiRs; sequences sharing very few mismatches with a miRNA and differing in length at the 5 and/or 3 ends) [47], the reads with more than two mismatches were considered as novel miRNAs in OSW. We propose to prepend the prefix 'asp' to the canonical miRNA name in order to associate it to OSW (e.g., asp-bantam, asp-miR-276a, asp-miR-317, etc.). IDs not present among the miRBase miRNAs were used for the novel miRNAs identified in this study, similarly preceded by the 'asp' prefix. Per each miRNA, the mature and star sequences were identified, as well as the indication of their 5p or 3p location in the precursor. Also, genomic coordinates of the precursors were determined in the B. tabaci genome. Normalized read count was expressed as reads per million (RPM) of reads mapping to the B. tabaci genome.
miRDeep2 calculates the number of raw reads per each predicted miRNA as the sum of the mature sequence and its isomiRs. In order to distinguish between read counts of mature miRNAs and isomiRs, we used the isomiRID tool [48]. Sequencing reads, mature and star sequences of miRNAs as well as their precursor sequences were submitted to isomiRID. Only isomiRs with up to one mismatch and differing in length up to ±2 and ±5 nt at the 5 and 3 ends, respectively, were retained. The output was manipulated in Microsoft Excel Professional 2016 for an arbitrary numbering of isomiRs. Also, the isomiR sequences were blasted against the miRBase database to find perfect matches with mature miRNAs. In order to assess the attitude of certain miRNA loci to generate more or less isomiRs, two indices were defined: the isomiR diversity index (IDI) was the number of different isomiRs relative to the normalized read count of a miRNA, and the isomiR abundance index (IAI) was the normalized read count of all isomiRs of a miRNA (excluding the mature sequence read count) relative to the normalized read count of the same miRNA.

Aleurocanthus spiniferus Was Found for the First Time on Ailanthus altissima in Italy
More than 100 puparia mounted on microscope slides were scrutinized. Based on taxonomic characters, the insects were identified as A. spiniferus because the number of marginal teeth ranged from 7 to 11, the crenulation sums 205-242 teeth, the ratio of longest spines were (1:9)-(1:17), the wide fringe represented 17.3%-30% of puparia width, cephalic eye spot were weakly defined, placed closer to the third sub-marginal spines rather than to the second ones, and microscopic papillae were situated between the sub-marginal spines [1,4]. Furthermore, the morphological identification was confirmed by a molecular study recently carried out on the same OSW population infesting the tree of heaven [3].
To the best of our knowledge, this is the first report of A. spiniferus infesting A. altissima. In fact, several studies have documented a number of pathogens and pests on A. altissima, but no mention of OSW has been reported [9][10][11][12][13]. This would suggest a host shift of OSW to A. altissima, but it cannot be excluded that the insect has infested this plant at extremely low population densities for a long time and it has not been observed or documented until now. Such a probable host shift merits a detailed investigation because it would add precious information about insect biology and ethology. It could be also studied whether insect endosymbionts, known as modulators of insect fitness [14,15], are involved in the colonization of this new plant host. Moreover, the infestation by a so aggressive insect pest on the tree of heaven could have implications in the control of this alien and invasive plant. In our inspections, the tree of heaven was not particularly threatened by OSW, even under heavy infestation levels. Therefore, a possible contribution by OSW for the control of the tree of heaven cannot be hypothesized hitherto.
3.2. "Candidatus Portiera", Wolbachia, and Rickettsia Were Identified in the OSW Microbiome DGGE profiles of the 16S rDNA showed three different dominant bands. All samples shared a band that could be ascribed to the Wolbachia reference, a known endosymbiont of Aleyrodidae [49] and other insect species [50]. Another band occurring in all samples had no correspondence to any employed positive controls. Furthermore, four out of ten samples showed a band corresponding to the Rickettsia reference, another endosymbiont of several insect species [51].
Sequencing and BLAST analysis confirmed that the DNA sequences obtained from the excised bands matched Rickettsia spp. and Wolbachia spp. accessions, while those obtained from the band with no correspondence to the used references matched accessions of "Candidatus Portiera" Thao and Baumann (Gammaproteobacteria) (here referred to as "Ca. Portiera").
By using deep sequencing of sRNAs we were able to infer a more variegated microbiome composition ( Figure 2). Most of the reads or contigs were unassigned to a Taxon. In fact, Kraken was able to classify only 0.28% reads, and Kaiju only 0.31%; among the 1289 contigs assembled with Velvet, only 242 were assigned to taxonomic units using a blastn search. On the other hand, 28 out of the 41 contigs assembled with VelvetOptimiser were classified by blastn.
Some microbial taxa were detected by different pipelines, while others only by one of them. Interestingly, the Aleyrodidae endosymbiont "Ca. Portiera" was identified by all the pipelines, and it was the most represented genus in three of them (except VelvetOptimiser). This bacterium is a well-known primary endosymbiont of Aleyrodidae species, including Aleurocanthus woglumi Ashby [49] and B. tabaci [14,52,53].
The Kraken tool detected considerable reads matching the Wolbachia genus. Also, the genera Rickettsia, Pseudomonas, Frateuria, and other taxa were identified by Kraken.
Kaiju identified the genera Methylobacterium Patt et al., Mycobacterium Lehmann and Neumann, and Streptomyces Waksman and Henrici among the bacteria, as well as Basidiomycota and Ascomycota (Bipolaris genus) among the fungi. Moreover, the genus Fragillariopsis (O'Meara) Hustedt (Bacillarophyceae) was also well represented. Besides Rickettsia, none of the latter microorganisms have been reported as Aleyrodidae endosymbionts. We did not find evidence for the presence of other secondary (facultative) endosymbionts of Aleyrodidae such as Hamiltonella Moran et al. [14,15,54,55] and another 14 genera [49]. Some species of the found guild of microorganisms might survive on the leaf surface, near or on the OWS; others seem frankly host-plant inhabitants that were taken by OSW during feeding and eventually concentrated in the filter chamber of the insect [56].
By using the contig-based pipelines (Velvet and VelvetOptimiser followed by blastn search), a number of sequences were assigned to Insecta, including Bemisia genus and other Aleyrodidae. No sequences were assigned to Aleurocanthus genus, probably because less represented than Bemisia or other Aleyrodidae in the databases. Remarkably, some contigs assembled by Velvet were assigned to the genus Eretmocerus Mercet (Hymenoptera), an Aleyrodidae parasitoid [57], and to the sub-order Cucujiformia (Coleoptera: Polyphaga), that encompasses plant-eating beetles and insect-predators [58].
Deep sequencing of insect sRNAs revealed that no plant-infecting viruses were present in the population of OSW hosted by the tree of heaven. When used with a virus database, Kraken detected three viral taxa: Mason-Pfizer monkey virus (M-PMV, Retroviridae; 2726 reads), Enterobacteriaceae phage phiX174 (Microviridae; 451 reads), and Herpesviridae (37 reads). Escherichia virus phiX174 was also detected by the Velvet + blastn pipeline (7 out of 242 contigs). Very likely, these hits were artifacts due to the forcing of search against viral sequences. Retroviridae (positive single-stranded RNA genome) and Herpesviridae (double-stranded DNA genome) infect vertebrates [59][60][61], but they have not been reported in insects, even after large screening surveys based on NGS [24][25][26]62]. Moreover, phiX174 (positive single-stranded DNA) is routinely used as a positive control in DNA sequencing on Illumina platforms [63]. Deep sequencing of insect sRNAs revealed that no plant-infecting viruses were present in the population of OSW hosted by the tree of heaven. When used with a virus database, Kraken detected three viral taxa: Mason-Pfizer monkey virus (M-PMV, Retroviridae; 2726 reads), Enterobacteriaceae phage phiX174 (Microviridae; 451 reads), and Herpesviridae (37 reads). Escherichia virus phiX174 was also detected by the Velvet + blastn pipeline (7 out of 242 contigs). Very likely, these hits were artifacts due to the forcing of search against viral sequences. Retroviridae (positive single-stranded RNA genome) and Herpesviridae (double-stranded DNA genome) infect vertebrates [59][60][61], but they have not been reported in insects, even after large screening surveys based on NGS [24][25][26]62]. Moreover, phiX174 (positive single-stranded DNA) is routinely used as a positive control in DNA sequencing on Illumina platforms [63]. The size of the circles is proportional to the number of reads or contigs placed in the taxon. Small RNA reads were subjected to Kraken and Kaiju or were assembled into contigs using Velvet or VelvetOptimiser prior to search sequence similarities against the BLAST nr database.

Features of Small RNAs
After quality filtering, we obtained 10,122,023 reads that were classified as tRNA, rRNA, miRNA, repetitive element, other Rfam families, others in B. tabaci genome, and unassigned RNA ( Figure 3B). About 29% of reads were tRNAs, and 41% remained unassigned, while the remaining reads were classified as rRNAs, miRNAs, repetitive elements, other Rfam elements, or other features in the B. tabaci genome. About 54% of reads (5,499,517) were mapped to the B. tabaci genome. We retained only reads with length ranging from 15 to 35 nt, and their length distribution is shown in Figure 3. Frequency peaks were observed for sRNAs of 21 and 28 nt, in agreement with previous works on B. tabaci B and Q biotypes or other insects [64][65][66]. Thirty-six percent reads between 19 and 22 nt mapped to miRNAs, while 74% of 28 nt reads were unassigned ( Figure 3B). They accounted for 8% of library ( Figure 4A), and 12% of the reads matching the B. tabaci genome ( Figure 3B; 5,499,517). Repetitive elements, accounting for 1% of library ( Figure 4A), were assigned to reads of a wide range Figure 2. Dendrograms of the microbiome of Aleurocanthus spiniferus infesting Ailanthus altissima based on the results obtained using four different tools for taxonomic assignment of small RNA reads. The size of the circles is proportional to the number of reads or contigs placed in the taxon. Small RNA reads were subjected to Kraken and Kaiju or were assembled into contigs using Velvet or VelvetOptimiser prior to search sequence similarities against the BLAST nr database.

Features of Small RNAs
After quality filtering, we obtained 10,122,023 reads that were classified as tRNA, rRNA, miRNA, repetitive element, other Rfam families, others in B. tabaci genome, and unassigned RNA ( Figure 3B). About 29% of reads were tRNAs, and 41% remained unassigned, while the remaining reads were classified as rRNAs, miRNAs, repetitive elements, other Rfam elements, or other features in the B. tabaci genome. About 54% of reads (5,499,517) were mapped to the B. tabaci genome. We retained only reads with length ranging from 15 to 35 nt, and their length distribution is shown in Figure 3. Frequency peaks were observed for sRNAs of 21 and 28 nt, in agreement with previous works on B. tabaci B and Q biotypes or other insects [64][65][66]. Thirty-six percent reads between 19 and 22 nt mapped to miRNAs, while 74% of 28 nt reads were unassigned ( Figure 3B). They accounted for 8% of library ( Figure 4A), and 12% of the reads matching the B. tabaci genome ( Figure 3B; 5,499,517). Repetitive elements, accounting for 1% of library ( Figure 4A), were assigned to reads of a wide range of length, but especially to the 17-18mers ( Figure 3B). The families L1-14_BDi, Gypsy-37_NVi-I, Gypsy-97_GM-I, Caulimovirus-23_ATr, and Gypsy-22_RO-I were the top-five represented transposable elements ( Figure 5).

MicroRNAs
Using the B. tabaci genome and miRBase database, the miRDeep2 tool identified 102 miRNAs in the OSW library. Sequence similarity analysis conducted with the blastn algorithm on these miRNAs against the miRBase database allowed us to classify the sequences into canonical (56) and novel miRNAs (46) (Tables 1 and 2, respectively; also see Tables S1 and S2 for more details). Reads sharing up to two mismatches with the database sequences were considered as canonical OSW miRNAs and the rest as novel miRNAs. Differently, a threshold of four mismatches was used in two previous works on B. tabaci, where more miRNAs than in our research were identified [64,65,67].
Canonical miRNAs were identified in 41 out of 19,751 scaffolds of the B. tabaci genome. MicroRNA precursors were predicted by miRDeep2, and 21 of them were highly similar to those available in miRBase (e.g., dme-mir-276a, dqu-mir-317, dme-mir-184, etc.). Several miRNAs were predicted to originate from diverse precursors, either with different nt sequences or with different positions in the genome. However, we cannot conclude if all these miRNAs actually have different precursors in OSW because we used the B. tabaci genome, and even the complete assembly of this genome into chromosomes is still pending. We found that asp-miR-iab-4 and asp-miR-iab-8 originated from the minus and plus DNA strand, respectively, of the same miRNA locus (Table 1), as already reported in D. melanogaster [68,69]. A similar scenario occurred for asp-miR-134420, which originated the same miRNA/miRNA* from both DNA strands. For most canonical miRNAs, we

MicroRNAs
Using the B. tabaci genome and miRBase database, the miRDeep2 tool identified 102 miRNAs in the OSW library. Sequence similarity analysis conducted with the blastn algorithm on these miRNAs against the miRBase database allowed us to classify the sequences into canonical (56) and novel miRNAs (46) (Tables 1 and 2, respectively; also see Tables S1 and S2 for more details). Reads sharing up to two mismatches with the database sequences were considered as canonical OSW miRNAs and the rest as novel miRNAs. Differently, a threshold of four mismatches was used in two previous works on B. tabaci, where more miRNAs than in our research were identified [64,65,67].    Canonical miRNAs were identified in 41 out of 19,751 scaffolds of the B. tabaci genome. MicroRNA precursors were predicted by miRDeep2, and 21 of them were highly similar to those available in miRBase (e.g., dme-mir-276a, dqu-mir-317, dme-mir-184, etc.). Several miRNAs were predicted to originate from diverse precursors, either with different nt sequences or with different positions in the genome. However, we cannot conclude if all these miRNAs actually have different precursors in OSW because we used the B. tabaci genome, and even the complete assembly of this genome into chromosomes is still pending. We found that asp-miR-iab-4 and asp-miR-iab-8 originated from the minus and plus DNA strand, respectively, of the same miRNA locus (Table 1), as already reported in D. melanogaster [68,69]. A similar scenario occurred for asp-miR-134420, which originated the same miRNA/miRNA* from both DNA strands. For most canonical miRNAs, we detected both mature and star sequences (distinguished based on their counts), supporting the reliability of those miRNAs in OSW according to the miRBase guidelines [45,46].
MicroRNA bantam promotes cell proliferation, inhibits apoptosis, maintains germline stem cells, limits scaling growth of dendrites, and regulates circadian rhythm clock [66,71,72]. MicroRNA miR-8 prevents neurodegeneration in the brain, is involved in the development of wings and eyes, promotes organismal growth in larvae, and controls the morphology of neuromuscular junctions [71]. MicroRNA miR-317 is critical for Drosophila survival, as it is a negative regulator of Drosophila Toll signaling immune response [73]. MicroRNA miR-750 is a cluster conserved in several hemimetabolan and holometabolan insects but lost in D. melanogaster [70]. It was the most abundant miRNA in silkworm infected by Bombyx mori cytoplasmic polyhedrosis virus (BmCPV) [67]. MicroRNA miR-276-5p fine-tunes the duration of the reproductive cycle in Anopheles [74] and promotes egg-hatching synchrony in L. migratoria [75]. Overall low expression of miRNAs, however, does not mean that they are of minor importance in the metabolism. MicroRNA mir-iab-4, for example, causes a dominant homeotic transformation of halteres to wings in D. melanogaster [76]. We found only 4.7 RPM of asp-miR-iab-4-5p, and miRBase reports 438 RPM as an average of 45 experiments for dme-miR-iab-4-5p. The abundance of some OSW miRNAs deviated from the general trend occurring in other insects. For example, let-7 had an intermediate abundance in OSW (745.2 RPM for the mature sequence) but was highly expressed in B. tabaci (1-2 × 10 5 ) and induced by TYLCCNV infection [64,65]. MicroRNA let-7-5p, along with miR-100-5p and miR-125-5p, emerged as an important factor in the metamorphic stage of B. germanica [77].
The tool miRDeep2 counts both the mature sequence and its isomiRs for each detected miRNA and miRNA star. In order to count separately the mature sequences and the isomiRs, we used the tool isomiRID and then extracted the counts from its output. Due to their abundance, isomiRs appeared particularly important to assess the expression level of a miRNA. For example, asp-bantam, which was the most abundant miRNA in the OSW library, accounted for 115,692.2 RPM (isomiRID count, Table S4), resulting from 33,695 RPM of mature miRNA and 79,930.3 RPM of isomiRs (606 different sequences). The isomiR diversity index (IDI), here defined as the number of different isomiR sequences relative to the mature miRNA count, decreased with the expression level of miRNAs, meaning that highly abundant miRNAs had a smaller number of different isomiRs ( Figure S1). However, such a correlation was not strong (R 2 = 0.42), so that miRNAs like asp-miR-2a-3p had 281.9 RPM of mature miRNAs and 181 different isomiRs (IDI = 0.64), while miRNAs like asp-miR-87a-3p had 240.57 RPM of mature miRNA and 31 isomiR sequences (IDI = 0.13). The isomiR abundance index (IAI), here defined as the total isomiR count relative to the mature miRNA count, did not depend on the miRNA abundance. In fact, miRNAs with similar abundance such as the two above-mentioned miRNAs asp-miR-2a-3p and asp-miR-87a-3p had a different IAI (1249.29 and 114.16 RPM of isomiRs, respectively, corresponding to IDIs of 4.43 and 0.47). In contrast, miRNAs with largely different expression levels such as asp-miR-317-3p (811.34 RPM of mature miRNA) and asp-miR-305-5p (99.1 RPM) had similar IDIs, i.e., 3.81 and 3.67, respectively.
By analyzing the isomiRs, we observed that A to C substitution was the most abundant occurrence (data not shown), which is not in line with Naqvi et al. [78], who found that most of the substitutions comprised A to G and C to T events.

Conclusions
We report the first documented stable infestation of A. altissima by OSW in Italy, which very likely originated from the neighboring citrus trees, the main host plants of OSW. A possible contribution by OSW for the control of the tree of heaven cannot be hypothesized hitherto because, in our observations, the plant tolerated well the OSW infestation.
By using DGGE and sRNA deep sequencing, we provide sound evidence on the presence of the primary endosymbiont "Ca. Portiera" and other endosymbionts such as Wolbachia spp. and Rickettsia spp. For the first time, we used sRNA-Seq to profile an insect microbiome and paved the way to compare this technique with the most widely used metagenomics approaches. Moreover, further investigations are needed to assess whether endosymbionts play a role in the interaction between OSW and the tree of heaven and whether the OSW microbiome can be harnessed to develop new management strategies against either the OSW or the tree of heaven.
Finally, we identified 56 canonical and 46 novel miRNAs in OSW, the latter group including at least four miRNAs detected with a high confidence level.
Supplementary Materials: The following are available online at http://www.mdpi.com/2075-4450/11/3/161/s1, Figure S1: Relationships between mature miRNAs and isomiRs of Aleurocanthus spiniferus. The isomiR diversity index is defined as the number of different isomiR sequences relative to the mature miRNA count (RPM). The isomiR abundance index is defined as the total isomiR count relative to the mature miRNA count (RPM). Table S1: Canonical miRNAs detected by miRDeep2 tool using Aleurocanthus spiniferus small RNA reads, Bemisia tabaci genome, and miRBase database. Table S2: Novel miRNAs (more than two mismatches shared with miRBase sequences) detected by miRDeep2 tool using Aleurocanthus spiniferus small RNA reads, Bemisia tabaci genome, and miRBase 22 database. Table S3: List and read count of mature miRNAs and isomiRs of Aleurocanthus spiniferus. Table S4: Overview of the isomiRs associated with the miRNAs of Aleurocanthus spiniferus.
Author Contributions: G.B. found the insect infestation, analyzed the data, and wrote the manuscript. F.P. and M.J. identified the insect species based on the morphological traits. G.B. and F.P. conceived the study. M.I.P. and F.G. collected the insects, handled them in the laboratory, and extracted the RNA. F.N. collected the insects and performed the PCR-DGGE assay. All the authors revised the manuscript. All authors have read and agreed to the published version of the manuscript.