Regulation of soldier caste differentiation by microRNAs in Formosan subterranean termite (Coptotermes formosanus Shiraki)

Methoprene bioassay

Termites (C. formosanus) were collected using a termite ground trap in Shangchong Fruit Tree Park, Guangzhou, China (Su & Scheffrahn, 1986; Chouvenc, Ban & Su, 2022). Only workers and soldiers were collected with the trap. The termites were then kept at room temperature and were used within two weeks of collection.

Soldiers were artificially induced using a JHA (methoprene) as follows. Filter paper (32 mm in diameter) was treated with acetone solution of methoprene and air dried for 15 min to evaporate the acetone, making the concentration of methoprene 1,000 ppm. The acetone-treated filter paper was used as a control. Two pieces of treated filter paper were placed in a plastic Petri dish (35 mm in diameter). A volume of 245 μL water was added to the two pieces of filter paper. Twenty workers were then put into the Petri dish. One colony of termites was used in this study. There were three replicates each in the methoprene-treated filter paper group (M group) and acetone-treated filter paper group (C group). After feeding on the filter paper for four days, the workers were collected and decapitated with a scalpel. The worker heads from the same Petri dish were placed into an Eppendorf tube with 100 μL Trizol reagent. The tube was then immersed in liquid nitrogen. All the samples were stored at −80 °C before use.

RNA extraction and sequencing

Total RNA was extracted using TRIzol reagent. The concentration and integrity of the RNA were measured using an Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA, USA). RNAs of 18–30 nt were enriched by polyacrylamide gel electrophoresis. Six libraries (three for the M group, three for the C group) were constructed by the NEBNext small RNA library prep set (New England Biolabs; E7300). Briefly, the 3’ and 5’ adapters were added to the RNAs. Then, the RNAs were reverse transcribed and amplified by PCR. The PCR products (140–160 bp in length) were recycled for sequencing. The libraries were evaluated by the Agilent 2100 and the ABI StepOnePlus Real-Time PCR System (Life Technologies, CA, United States). Finally, the libraries were sequenced by Illumina Hiseq2000 (Illumina, San Diego, CA, USA).

Identification of miRNAs and miRNA target prediction

The sequence of the adapters and the low-quality bases were filtered out to acquire clean tags. The following reads were filtered out: reads containing more than one low-quality (Q-value ≤ 20) or unknown nucleotides; reads without 5’ or 3’ adapters; reads without a nucleotide between the 5’ and 3’ adapters; reads containing polyA; and reads shorter than 18 nt (not including adapters). The clean reads were then aligned with small RNAs in GeneBank and Rfam using Blastall 2.2.25 (blastn). The reads of the non-coding RNAs (ncRNAs), including rRNAs, scRNAs, snoRNAs, snRNAs, and tRNAs, were removed from the clean tags. The miRNAs of C. formosanus were not included in the miRBases, so the remaining tags were searched against the miRBases database to identify known miRNAs by aligning with the miRNAs of other species using Bowtie (version 1.1.2) (Langmead & Salzberg, 2012). For the known miRNA of C. formosanus, x denoted that the miRNA was aligned with a 5’ mature miRNA in the miRBases database; y denoted that the miRNA was aligned with a 3’ mature miRNA in the miRBases database. Possible novel miRNAs were identified according to their position on the unigenes and hairpin structures predicted by Mireap_v0.2 using the default parameters for animals. The reads of the known miRNAs and novel miRNAs were then removed from the clean tags, and the remaining tags were aligned with the RNA-Seq data of C. formosanus to remove the reads from mRNA degradation. The RNA-Seq data were obtained by sequencing the transcriptomes of C. formosanus workers feeding on methoprene-treated filter paper and workers feeding on acetone-treated filter paper for four days (Du et al., 2020). After identification of ncRNAs, known miRNAs, novel miRNAs, and reads from mRNA degradation, the rest of the tags were recorded as unannotated.

The transcriptomes of C. formosanus were used for miRNA target prediction (Du et al., 2020). The target genes of miRNAs were predicted by RNAhybrid (v2.1.2) + svm_light (v6.01), Miranda (v3.3a), and TargetScan (Version:7.0). The intersection of the target genes identified by all three of these methods were identified as the miRNA target genes.

Analysis of sample relationship and differentially expressed miRNAs (DEmiRs)

The expression levels of miRNAs were normalized to transcripts per million (TPM). The normalized expression levels of miRNAs were subjected to a principal component analysis (PCA) using the R package gmodels (R Core Team, 2013). The software edgeR was used to identify the differentially expressed miRNAs between the C group and the M group with the criteria |log₂(fold change)| > 1 and P value < 0.05 (Robinson, McCarthy & Smyth, 2010). The statistical power analysis was performed using RNASeqPower.

Validation of miRNA sequencing results

The expression levels of the miRNAs were quantified using the methods outlined by Chen et al. (2005). RNA samples were reverse transcribed using stem-loop RT primers and the cDNAs were analyzed by dye-based qPCR. The stem-loop RT primers for the target genes and the reference gene were used to reverse transcribe the RNA samples with M-MuLV (NEB, Ipswich, Massachusetts, US). The PowerUp SYBR Green Master Mix (Thermo Fisher Scientific, Waltham, MA, USA) was used to conduct the qPCR with the miRNA specific forward primer and a universal reverse primer (sequence: GTGCAGGGTCCGAGGT). Four miRNAs were selected from the DEmiRs, i.e., let-7-y, miR-1175-y, miR-181-x, and miR-305-y. Novel-m0649-y was used as the reference gene based on the screening results of the most stable reference gene (Du et al., 2023). The sequences of the stem-loop RT primers and the forward primers are summarized in Table S1. The data were analyzed using the 2^−ΔΔCt method (Livak & Schmittgen, 2001). The 2^−ΔΔCt data were log₂ transformed and analyzed using a t-test to compare the expression value between the C group and the M group.

Determination of DEmiR-regulated genes

The differentially expressed genes (DEGs) were obtained from the transcriptome data of C. formosanus (Du et al., 2020). Because miRNAs mediate the repression of mRNA expression post-transcriptionally, the DEmiRs with the potential to control soldier caste differentiation should have an inverse expression relationship with their targeting DEGs. For visualization of the DEmiR-DEG regulatory network, a pair-wise Spearman correlation was conducted for DEmiR-DEG targeting pairs using R. The DEmiR-DEG targeting pairs with Spearman’s rank correlation coefficient <−0.6 and P < 0.05 were selected to construct the DEmiR-DEG regulatory network with Cytoscape. For clarity, only the top 18 most expressed DEmiRs were shown from a total of 116 DEmiRs.

Enrichment analysis of the DEmiR-targeted DEGs

Gene Ontology (GO) systematically describes the features of genes and gene products with a standardized vocabulary in three categories: cellular components, molecular functions, and biological processes. The Kyoto Encyclopedia of Genes and Genomes (KEGG) is a database incorporating molecular-level data to assist in the exploration of high-level functions and utilities of biological systems. KEGG pathway establishes a relatively comprehensive group of pathway maps about interaction, reaction, and relation networks at the molecular level. In order to evaluate the functions of the DEmiRs in soldier caste differentiation and identify the significant GO terms and KEGG pathways that are regulated by miRNAs, the DEmiR-regulated DEGs were subjected to the GO and KEGG pathway enrichment analysis using Fisher’s exact test. False discovery rate was used to adjust the P value.

Sequencing quality

Six libraries were constructed, three for each treatment. The sequencing quality of the libraries is shown in Table S2. There were 74,791,183 clean reads after filtering out the invalid sequences. After removing the low-quality reads and contaminants, including adapters, inserts, polyA reads, reads smaller than 18 nt, and low cutoff reads, a total number of 67,959,005 clean tags were acquired. There were 33,598,522 and 34,360,483 clean tags for the three libraries in the C group and the M group, respectively (Tables 1 and S2). All the samples had more than 91% clean tags. The length distribution of the tag is a characteristic animal small RNA, with the highest percentage at 22 nt. The percentage of the 22 nt tag was 33.93%, 33.91%, 34.41%, 34.29%, 33.37%, 34.97% for C-1, C-2, C-3, M-1, M-2, M-3, respectively (Fig. S1).

Identification of miRNAs and target prediction

Four ncRNAs, including rRNAs, snRNAs, snoRNAs, and tRNAs, were identified by searching the Rfam ncRNA database, while only two ncRNAs (rRNAs and tRNAs) were identified using the Genbank ncRNA database. However, more counts of rRNAs and tRNAs were identified using Genbank than when using Rfam (Table S3). The percentage of known miRNAs was the highest among all the tags in either the M or the C group, and its percentage in the C group was similar to that in the M group (76.52 ± 0.47% vs 77.37 ± 1.93%; Table 1). However, the reads annotated as novel miRNAs only accounted for a small portion of the total reads (0.39 ± 0.01% in the C group, 0.37 ± 0.00% in the M group). There were 6.82 ± 0.19% counts which were identified as transcriptome data in the C group, and 6.24 ± 0.08% in the M group. A total of 1,324 miRNAs were identified, with 589 known miRNAs and 735 novel miRNAs (Table 2). In addition, 21,281 genes were predicted as the targets of the miRNAs. The miRNAs and their computationally-predicted-targets are summarized in Table S4.

Sample relationship between treatments

The PCA analysis showed high variability among samples (Fig. 1). The first component only accounted for 35.6% of the variance, while the second component accounted for 20.2% of the variance. According to the 3-dimensional figure of the PCA analysis, the third component accounted for 17.3% of the variance. The three components together accounted for 73.1% of the total variance.

DEmiRs between treatments

There were 116 differentially expressed miRNAs between the C group and the M group: 16 up-regulated miRNAs and 100 down-regulated miRNAs in the M group, using the C group as the reference. Some of the DEmiRs were expressed at relatively low expression levels. The DEmiRs, their sequences, and their expression levels are shown in Table S5. The statistical power of this experimental design was 0.41 (n = 3).

Validation of miRNA sequencing results

The qPCR data were only partially consistent with the miRNome data. There were significant differences between the C group and the M group for let-7-y and miR-305-y (Figs. 2A and 2B). According to the qPCR results, there was no significant difference for miR-1175-y and miR-181-x between the C group and the M group (Figs. 2C and 2D), but the miRNome sequencing data showed different expression levels of miR-1175-y and miR-181-x between the C group and the M group.

DEmiR-DEG pairs and functional analysis of the DEmiR-regulated DEGs

A total of 2,547 DEGs were found between the C group and the M group after treating workers with methoprene for four days, including 1,480 up-regulated genes, and 1,067 down-regulated genes in the M group, using the C group as the reference (Du et al., 2020). In combination with the results of the inverse relationship and the targeting relationship, a total of 4,433 DEmiR-DEG pairs were obtained (Table S6). The DEGs in the DEmiR-DEG pairs were considered to be the genes regulated by the DEmiRs. Among the top 18 most expressed DEmiRs, miR-148-y targeted the most DEGs that were down-regulated after workers were treated with methoprene (Fig. 3A), while let-7-y targeted the largest number of genes that were up-regulated after methoprene treatment (Fig. 3B). The KEGG pathway classification and GO annotation of the DEmiR-targeted DEGs are shown in Table S6. The DEmiRs were involved in the regulation of genes in metabolic pathways, the TGFβ signaling pathway, hedgehog signaling pathway, and FoxO signaling pathway (Table S6).

Although several GO items in the cellular component, molecular function, and biological process categories were enriched according to the P values, no GO term was recognized as significant after adjusting the P values (Fig. 4). The KEGG enrichment analysis showed that there were two significantly enriched pathways: ribosome biogenesis in eukaryotes (Q < 0.009) and the circadian rhythm-fly pathway (Q < 0.036; Fig. 5). The top 10 enriched pathways were mainly related to biogenesis and metabolism, with the circadian rhythm-fly pathway being the only enriched KEGG pathway related to signal transduction.