Identiﬁcation and Characterization of Development-Related microRNAs in the Red Flour Beetle, Tribolium castaneum

: MicroRNAs (miRNAs) play important roles in insect growth and development, but they were poorly studied in insects. In this study, a total of 883 miRNAs were detected from the early embryo (EE), late larva (LL), early pupa (EP), late pupa (LP), and early adult (EA) of Tribolium castaneum by microarray assay. Further analysis identiﬁed 179 differentially expressed unique miRNAs (DEmiRNAs) during these developmental stages. Of the DEmiRNAs, 102 DEmiR-NAs exhibited stage-speciﬁc expression patterns during development, including 53 speciﬁcally highly expressed miRNAs and 20 lowly expressed miRNAs in EE, 19 highly expressed miRNAs in LL, 5 weakly expressed miRNAs in EP, and 5 abundantly expressed miRNAs in EA. These miRNAs were predicted to target 747, 265, 472, 234, and 121 genes, respectively. GO enrichment analysis indicates that the targets were enriched by protein phosphorylation, calcium ion binding, sequence-speciﬁc DNA binding transcription factor activity, and cytoplasm. An RNA interference-mediated knockdown of the DEmiRNAs tca-miR-6-3p, tca-miR-9a-3p, tca-miR-9d-3p, tca-miR-11-3p, and tca-miR-13a-3p led to defects in metamorphosis and wing development of T. castaneum . This study has completed the identiﬁcation and characterization of development-related miRNAs in T. castaneum , and will enable us to investigate their roles in the growth and development of insect.


Introduction
microRNAs (miRNAs) are a type of endogenous non-coding small RNAs that can regulate post-transcriptional gene expression through complementary base-pairing with the messenger RNA (mRNA) of the target gene to form RNA-induced silencing complex [1]. Usually, primary miRNA is transcribed by RNA polymerase-II and then cleaved by the Drosha-DGCR8 complex to form a smaller structure called pre-miRNA [2]. This structure is transported to the cytoplasm by Exportin-5 and then enzymatically cleaved by Dicer to form mature miRNA [3,4]. miRNAs are associated with many biological processes, such as development, metabolism, apoptosis, autophagy, immunity and virus transmission [5][6][7][8][9][10].

Identification of DEmiRNAs between Two Developmental Stages
Since miRNAs from embryos, larvae, pupae, and adults have been obtained, it is necessary to identify DEmiRNAs between the two developmental stages. By microarray

Identification of DEmiRNAs between Two Developmental Stages
Since miRNAs from embryos, larvae, pupae, and adults have been obtained, it is necessary to identify DEmiRNAs between the two developmental stages. By microarray analysis, we found 142 DEmiRNAs (97 upregulated miRNAs and 45 downregulated miRNAs) between LL and EE, 140 DEmiRNAs (47 and 93) between EP and LL, 76 DEmiRNAs (50 and 26) between LP and EP, and 11 DEmiRNAs (7 and 4) between EA and LP ( Figure S1). From the volcano plots of sequences with strong signal (≥500), 73, 84, 47, and 4 DEmiRNA sequences were found in pairs of LL-VS-EE, EP-VS-LL, LP-VS-EP, and EA-VS-LP, respectively (p < 0.05). Of these, 54 miRNA sequences were upregulated whereas 19 were downregulated during EE to LL development. Thirty-eight miRNAs were upregulated while 46 were downregulated when LL molt into EP. Thirty-three miRNAs were upregulated whereas 14 were downregulated during EP to LP development. One miRNA was upregulated and 3 miRNAs were downregulated when LP entered to EA (Figure 3).

Analysis of Stage-Specific DEmiRNAs
In order to explore the role of miRNAs in specific developmental stage, we next conducted the identification of stage-specific miRNAs. According to the expression patterns, 102 DEmiRNAs with the signal strength ≥500 were classified into five clusters ( Figure 4A-E). The cluster A includes 53 embryo-specific highly expressed miRNAs, the cluster B comprises 20 embryo-specific lowly expressed miRNAs, the cluster C contains 19 miRNAs specifically highly expressed in LL, the cluster D consists of 5 miRNAs with specifically weak transcription in EP, and the cluster E has 5 miRNAs specifically highly expressed in EA (Table 1).

Analysis of Stage-Specific DEmiRNAs
In order to explore the role of miRNAs in specific developmental stage, we next conducted the identification of stage-specific miRNAs. According to the expression patterns, 102 DEmiRNAs with the signal strength ≥ 500 were classified into five clusters ( Figure 4A-E). The cluster A includes 53 embryo-specific highly expressed miRNAs, the cluster B comprises 20 embryo-specific lowly expressed miRNAs, the cluster C contains 19 miRNAs specifically highly expressed in LL, the cluster D consists of 5 miRNAs with specifically weak transcription in EP, and the cluster E has 5 miRNAs specifically highly expressed in EA (Table 1).

Target Prediction of DEmiRNAs
MiRNAs usually bind to targets to participate in multiple biological processes, and we subsequently predicted the targets of these DEmiRNAs by miRanda, PicTar, and Tar-getScan, respectively. A total of 2025 target genes with 18,734 target sites, 2953 target genes with 30,009 target sites, and 2205 target genes with 24,204 target sites were screened from miRanda, PicTar, and TargetScan, respectively. There are 68.4% common spots between miRanda and PicTar, 63% common spots between miRanda and TargetScan, and 74.4% common spots between PicTar and TargetScan. All three softwares shared 63% common spots. Totally 1863 common target genes were used for further analysis ( Figure 5). For the DEmiRNAs, 747, 265, 472, 234, and 121 target genes were obtained from EE-specific highly expressed miRNAs (cluster A) and lowly expressed miRNAs (cluster B), LL-specific highly expressed miRNAs (cluster C), EP-specific lowly expressed miRNAs (cluster D), and EA-specific highly expressed miRNAs (cluster E), respectively (Table 1).

DEmiRNAs Are Required for Metamorphosis and Wing Development
Among these miRNAs, we selected five miRNAs (three from cluster A, and one each of cluster C and E) that have been revealed by GO analysis to play critical roles in development of the insect to verify the relationship between these miRNAs and development.

Discussion
In this study, we identified DEmiRNAs among five developmental stages by microarray, predicted their targets, functionally annotated them by GO analysis of their targets, and investigated the developmental roles of five DEmiRNAs in T. castaneum.
In the past years, a number of miRNAs have been identified from T. castaneum [15,18]. By comparison with other developmental stages, we found 53 specially upregulated miR-NAs (cluster A) and 20 downregulated miRNAs (cluster B) in the early embryonic (EE) stage of T. castaneum (Table 1, Figure 4A,B). Among the upregulated miRNAs, miR-7 exerts a role in the switch from the endocycle to gene amplification through its regulation of Tramtrack69 in Drosophila follicle cells [44], miR-9 specifically controls the expression of mesodermal genes to canalize myotendinous junction morphogenesis during embryogenesis [45], and miR-309 targets SIX4 and controls ovarian development in A. aegypti [46]. Of the downregulated miRNAs, let-7 could bind to the Kr-h1 coding sequence to block oocyte maturation in Locusta migratoria [47], miR-1-3p is an early embryonic male sexdetermining factor in the Oriental fruit fly B. dorsalis [48], miR14 regulates egg-laying by targeting EcR in A. mellifera [49], loss of miR-184 in Drosophila led to multiple severe defects during oogenesis and early embryogenesis, culminating in the complete loss of egg production [50], and miR-276 promotes egg-hatching synchrony by upregulating brahma in locusts [51]. These studies indicate that specifically high or low expression of miRNAs in EE is associated with embryogenesis of T. castaneum.
This study discovered 19 miRNAs specifically highly expressed in the LL of T. castaneum (cluster C) (Table 1, Figure 4C). Of these miRNAs, an overexpression of miR-8 in the whole body at the end of Drosophila larval development inhibited ecdysone biosynthesis and caused defects in metamorphosis and survival [52], and its overexpression in the corpus allatum increased cell size of the gland and expression of JHAMT, leading to pupal lethality [37]. Specific inhibition of miR-8-3p in T. castaneum late larvae induced metamorphosis defects in the development of wings, eyes, legs and embryos [28]. The ablation of B. mori miR-34 by transgenic CRISPR/Cas9 modulated ecdysone signaling and resulted in a severe developmental delay during the larval stage [53], and miR-281 regulates the expression of EcR isoform B in B. mori [54]. The injection of tca-miR-6-3p antagomirs caused defects in ecdysis during pupation and eclosion and wing development ( Figure 8). Therefore, these miRNAs were believed to have critical roles in metamorphosis of T. castaneum.
The cluster D has 5 miRNAs specifically weakly expressed in EP (Table 1, Figure 4D), and they were predicted to target 234 protein-coding genes ( Table 1). The GO enrichment analysis showed that these DEmiRNAs mainly participate in the cytoplasm (32), plasma membrane (20), signal transduction (15), and sequence-specific DNA binding transcription factor activity (14) ( Figure 7B). In the cluster, miR-9a targeting senseless ensures the precise specification of sensory organ precursors in Drosophila [55], and it also target the Drosophila transcriptional regulator LIM-only to prevent apoptosis during wing development [56]; miR-8 null flies are smaller in size and defective in insulin signaling in fat body [57], gain and loss of miR-8 provoked developmental defects reminiscent of impaired Notch signaling [58], and it also acts as a negative regulator of Wnt signaling to affect eye and wing development in D. melanogaster [59]; miR-252 has vital roles in growth, development and cell cycle of insect [9,60,61]; miR-3017 contributes to larval to pupal to adult metamorphosis by targeting sarco/endoplasmic reticulum Ca 2+ ATPase in T. castaneum [62]. Mutation of the Drosophila miR-6 cluster, which contains three miR-6 genes, showed a modest reduction in viability to adulthood [63]. A knockdown of tca-miR-6-3p also led to defects during the transition of larva to pupa to adult (Figure 8). These results suggest that the cluster D miRNAs have multiple functions in organ growth, development, and cell cycle.
The cluster E contains 5 EA-specific highly expressed miRNAs (Table 1, Figure 4E). Among these DEmiRNAs, the dietary of miR-316 antagomirs led to significantly higher mortality (>70%) and a lower proportion of winged morphs [64], miR-9 functions in muscle contraction [65], dendrite growth of sensory neurons [66], and wing development ( Figure 8) [67] of insect, and miR-375 has important roles in forming short-term memory [68], circadian rhythms and sleep [69] of D. melanogaster. Therefore, these miRNAs participate in many aspects of adult development in insects.
In conclusion, we identified 179 DEmiRNAs among five developmental stages of T. castaneum, and found 102 DEmiRNAs specifically expressed in EE, LL, EP, and EA. GO analysis indicated these miRNAs played vital roles in development. A knockdown of five DEmiRNAs caused defects in metamorphosis and wing development. This study completed the identification of development-related miRNAs in T. castaneum.

Insect Strains
The Georgia-1 (GA-1) strain of T. castaneum was reared at 30 • C in 5% yeast flour [70]. They were used in all experiments of this study.

RNA Isolation and cDNA Synthesis
The total RNA was extracted from early embryos (1-day embryos), late larvae (20-day larvae), early pupae (1-day pupae), late pupae (6-day pupae), and early adults (1-day adults) using a Total RNA Purification Kit (LC Sciences, Houston, TX, USA) according to the manufacturer's protocol, respectively. The quantity and purity of the total RNA were analyzed with the Bioanalyzer 2100 by RNA 6000 Nano LabChip Kit (Agilent, Houston, TX, USA) with RIN number >7.0. In addition, 500 ng of total RNA was converted to cDNA using HiScript ® III RT SuperMix for qPCR (+gDNA wiper) (Vazyme, Nanjing, China).

Microarray Assay
A poly (A) tail was added to the 3 end of 4 µg total RNAs of each sample by the poly (A) polymerase. An oligonucleotide tag was then ligated to the poly (A) tail for subsequent fluorescent dye staining. Hybridization was performed overnight on a µParaflo microfluidic chip by a micro-circulation pump (Atactic Technologies, Texas, TX, USA) [71]. On the microfluidic chip, each detection probe consisted of a chemically modified nucleotide coding segment that was complementary to the target miRNA and a spacer segment of polyethylene glycol to extend the coding segment away from the substrate. The detection probes were made in situ synthesis using photogenerated reagent (PGR) chemistry. The hybridization melting temperatures were balanced by chemical modifications of the detection probes. Hybridization used 100 L 6 × SSPF buffer (0.90 M NaCl, 60 mM Na 2 HPO 4 , 6 mM EDTA, pH 6.8) containing 25% formamide at 34 • C. After RNA hybridization, tag-conjugating Cy3 dye was circulated through the microfluidic chip for dye staining. Fluorescence images were collected using a laser scanner (GenePix 4000B, Molecular Device, California, CA, USA) and digitized using Array-Pro image analysis software (Media Cybernetics, Maryland, MD, USA). Data were analysed by first subtracting the background and then normalizing the signals using a LOWESS filter (Locally-weighted Regression, California, CA, USA) [72]. Each microarray contained 16 groups of 5sRNA probes as control, and these controls of signal would >10 to insured small RNAs were existed in the detected samples.

Identification of Differentially Expressed miRNAs
Fluorescence images were collected using a laser scanner (GenePix 4000B, Molecular Device) after RNA hybridization. The "Array-Pro Analyzer" (Media Cybernetics) was used for image digitization and background subtraction [71]. Normalization was carried out using the locally weighted scatterplot smoothing (LOWESS) method on the subtracted background data [72]. Differentially expressed miRNAs (DEmiRNAs) between developmental stages were filtered using ANOVA and t-tests [73].

Target Prediction for DEmiRNAs
Four different algorithms, including NCBI blast, PicTar (http://pictar.mdc-berlin.de/, accessed on 20 May 2022), miRanda (http://mirdb.org/miRDB/, accessed on 20 May 2022), and TargetScan (http://www.targetscan.org/vert_61/, 20 May 2022), were used to predict the potential targets of DEmiRNAs. PicTar calculates the free energy of miRNA, which represents the binding capability to the target gene [74]. miRanda analyses the free energy of miRNA binding capability and the binding site between miRNA and target gene [75]. TargetScan searches for complementary matches (called the "seed match") of miRNA 2-8 bases (numbered from the 5 end) with the 3 UTR sequence of mRNA [76]. Only miRNA-target interactions predicted by all three algorithms were used for further analysis.

Enrichment Analyses of Predicted DEmiRNA Targets
We used Blast2GO [77] to annotate the gene ontology (GO) term for DEmiRNA targets. A GO term with false discovery rate (FDR) ≤ 0.01 was judged to be a significantly enriched one.

Quantitative Real-Time PCR
Total RNAs were separately collected from samples and then converted to cDNAs for each miRNA by RT-primer (Table S1) using HiScript ® III RT SuperMix for qPCR (+gDNA wiper) (Vazyme, Nanjing, China). The quantitative real time-PCR (qRT-PCR) was performed using the Hieff q-PCR Green Master Mix (Yeasen, Shanghai, China) to examine the expression of miRNAs. The relative levels of miRNAs were normalized to the level of U6 snRNA using the ∆∆ CT method [78]. The results were analyzed from three independent samples. The primer sequences for these miRNAs are listed in Table S1.

Statistical Analysis
Data were statistically analyzed by Student's t-test and analysis of variance (ANOVA) using the graphic software Prism (Graph Pad Software, v8.2.1, San Diego, CA, USA). All data are presented as the mean ± standard error of the mean (SEM).