An Aedes albopictus laboratory line was established from wild mosquitoes which were collected from the region of Thessaly, as described elsewhere (Ioannou et al., 2021), and were reared in the insectary facility of the Department of Biochemistry & Biotechnology at the University of Thessaly. Adult mosquitoes were reared at 26 ± 1°C with a relative humidity of 60−70%, under a 14 h:10 h light/dark photoperiod. They were fed with 10% sucrose solution, while female mosquitoes were blood-fed (BF) from a human arm to initiate their gonotrophic cycles.

Ovaries and other tissues were dissected under the microscope and their total RNA was extracted using Extrazol (BLIRT S.A., Gdańsk, PL), according to the manufacturer’s instructions. The integrity of the RNA was assessed by agarose gel electrophoresis. Total RNA was treated with DNase I (Thermo Fisher Scientific, Waltham, MA, United States) and 1 μg of RNA was reverse transcribed to cDNA by using oligo-dT primers and MMLV-RT (Invitrogen, Waltham, MA, United States), according to the manufacturer’s instructions. Each biological replicate corresponds to tissues collected from individual mosquitoes. We preferred to study tissues from individual mosquitoes, rather than pooling them together, in order to be able to assess within population variability. All qPCR assays were performed in CFX96 Real-Time Thermal Cycler (Bio-Rad, Hercules, CA, United States). All amplifications were performed with two technical replicates and the relative gene expression was analyzed by using the 2^−ΔΔCt method (Livak) through the CFX Manager™ software. Specific primers to amplify the genes identified by the transcriptomic analysis were designed using PrimerQuest™Tool. Their target specificity was verified through Primer-BLAST (Ye et al., 2012) against the RefSeq mRNA database of Aedes albopictus. Primers that lacked any sequence homology to other transcripts were selected. Two endogenous ribosomal housekeeping genes (RpL32, RpS17) were used for the normalization of the results (Dzaki and Azzam, 2018). The average expression stability value (M-value) of the reference genes for each biological sample was determined and samples that exhibited M-value <0.5 and Coefficient of Variance <0.25 were accepted for further analysis. Primers for qPCR are shown in Supplementary Table S2.

Target templates for in vitro transcription were generated using gene specific primers with the respective recognition site for T7 RNA polymerase designed by the eRNAi web platform (Horn and Boutros, 2010). Thermodynamic parameters of the primers were tested through the web platform OligoAnalyzer (https://eu.idtdna.com/calc/analyzer). We selected primer sets that exhibited a minimum ΔG value of −9 kcal/mol for homodimer and heterodimer formation, while their hairpin structures did not exceed the primer Tm-10°C. The target specificity of the primers was assessed by Primer-BLAST (Ye et al., 2012) against the RefSeq mRNA database of Aedes albopictus. The off-target effect of dsRNAs (T7 amplicons) that were introduced for silencing was evaluated by examining for potential homologies of their siRNAs with other transcripts of Aedes albopictus. We obtained all possible 21-mers (i.e., sliding window with width = 21 and step = 1) and subjected them to remote blastn analysis (argument -task “blastn-short” for short input queries) against RefSeq mRNA database of Aedes albopictus. We proceeded with dsRNAs/siRNAs that did not contain 21-mers matching any other transcript of Aedes albopictus, other than their lncRNA target. dsRNA was synthesized by using the MEGAscript RNAi kit (Ambion, Austin, TX) from a dsDNA template by incubating at 37°C for 16 h. dsDNA was produced by PCR with cDNA template and gene-specific primers with the T7 RNA polymerase promoter (TAA​TAC​GAC​TCA​CTA​TAG​GG) attached to their 5’ end. dsRNA was purified with standard phenol/chloroform protocol, its integrity was measured by agarose gel electrophoresis and its concentration was quantified by the Q3000 Spectrophotometer (Quawell, San Jose, CA). Green fluorescent protein (GFP) dsRNA was used as a control. For silencing with each one of the Norma genes, inseminated Ae. albopictus adult females were used. The mosquitoes were collected 5 days after eclosion and were anesthetized on CO₂. Then, 64.4 nl of dsRNA solution (5,000 ng/μl) was injected into their thorax using the Nanojet II microinjector (Drummond, Birmingham, AL), under a Leica™ stereoscope. At least 25 individual females (i.e., biological replicates) were used for silencing with each Norma gene. Two controls were used for comparison: one was non-injected females (untreated) and the other was GFP-dsRNA injected females (anti-GFP). Both control populations were reared under the same conditions with the Norma-dsRNA injected ones. Primers for dsRNA production are shown in Supplementary Table S2.

Five-day-old mosquitoes injected with either GFP-dsRNA or Norma3-dsRNA were transferred to Bugdorm™ 17.5 × 17.5 × 17.5 cm cages (MegaView Science Co., Talchung, Taiwan) immediately after injection, where they were fed with 10% sucrose solution for 24 h to recover. An additional sample of same age, non-injected mosquitoes (untreated) was also maintained. The mosquitoes were left to starve for 12 h and subsequently they were blood-fed (36 h post-injection). All blood meals took place in the same timeframe, between morning and noon, to avoid the influence of circadian variability among the replicates. Fully engorged, blood-fed mosquitoes were reared in the presence of 10% sucrose.

Raw FASTQ files from a detailed, publicly available, developmental transcriptomic dataset by Gamez and colleagues were retrieved from the Sequence Read Archive (SRA) study entry SRP219966 (Gamez et al., 2020), using SRA toolkit. Quality assessment, identification of adapters and overrepresented contaminants was performed by employing FASTQC (www.bioinformatics.babraham.ac.uk/projects/fastqc/) and the EMBL-EBI’s Minion application (http://www.dev.ebi.ac.uk/enright-dev/kraken/reaper/src/reaper-latest/doc/minion.html). Quality trimming and adapter trimming of reads was performed with cutadapt. Pre-processed FASTQ files were mapped to the GCF_001876365.2_canu_80X_arrow2.2 Ae. albopictus genome assembly utilizing STAR splice-aware aligner (v2.7.9a) (Dobin et al., 2013) using the ENCODE standard options provided in the STAR manual. The respective transcript annotation file, derived from the NCBI Eukaryotic Genome Annotation Pipeline (Gnomon gene predictions), was also provided (--sjdbGTFfile argument) to effectively account for known splice junctions, while the genome index was accordingly parameterized to the read length (--sjdbOverhang 49). Transcript-level expression was calculated with RSEM v1.3.1. (Li and Dewey, 2011).

Accession numbers with the “XR_” prefix, annotated as non-coding RNA (ncRNA) were retained from the GCF_001876365.2 gene models. In order to obtain estimates of the transcripts’ coding potential, CPC2 (Kang et al., 2017) and FEELnc (Wucher et al., 2017) tools were employed. Both tools were executed at default settings (in the absence of any known lncRNAs, FEELnc “--mode = shuffle” method was chosen and the “XM_” (mRNA) transcripts were also provided as input for training). Transcripts were queried locally with BLASTn v2.13.0 (Altschul et al., 1990), at default settings, against Reference RNA Sequences (RefSeq_rna) from all Hexapoda species (Taxonomy ID: 6960). Each transcript was annotated with respect to the number of species it presented hits to, besides Ae. albopictus. FEELnc classifier was also used to characterize all “XR_” transcripts as genic and intergenic, according to their genomic localization and strand, respective to other transcripts. Lastly, RSEM-derived counts from all analyzed datasets in Section 2.6 were TMM-normalized (Robinson and Oshlack, 2010), log₂-transformed and utilized to obtain tau (τ) tissue specificity indices (Yanai et al., 2005) (https://github.com/roonysgalbi/tispec). For a transcript x with expression x _i across n tissues/contexts, tau is obtained using its expression normalized by its maximal expression as follows: t a u = ∑ 1 n ( 1 − x ^ i ) n − 1   ;   x ^ i =   x i max ( x i )   1 ≤ i ≤ n

After deriving the global tau for a transcript, its per-context specificity fraction was calculated by multiplying tau with the maximal-normalized expression metric x ^ i . “XR_” transcripts were annotated with their context (i.e., developmental stage, tissue and time-point) tau specificity index, as well as per-context fraction. Heatmaps were created using R package pheatmap, using Euclidean distance as a metric and complete-linkage hierarchical clustering where applicable. Annotation results are provided unfiltered in Supplementary Data S1.

For correlation analysis, samples where Norma3 presented non-zero expression (Transcripts-Per-Million, TPM units) were selected and transcripts presenting zero TPM across all samples were filtered out. Pearson correlation testing and False Discovery Rate adjustment of p-values were conducted in R. Clustering analysis was performed by providing TPM values to DPGP (McDowell et al., 2018), which applies Dirichlet Process to nonparametrically determine the optimal number of expression trajectory clusters and Gaussian Process to model the trajectories of expression through time. Following DPGP recommendations, a limited list of transcripts was subjected to clustering analysis instead of the whole transcriptome. For that purpose, read counts from post-blood meal (PBM) ovary samples (12, 24, 36, 48, 60, and 72 h) were imported in R. In the absence of replicates, 12–24 h, 36–48 h and 60–72 h time points were grouped into “early”, “middle” and “late” groups respectively. EBSeq (v1.30.0) (Leng et al., 2013) was utilized to assign transcripts to patterns of differential expression, at a relaxed posterior probability ≥75%. DPGP was executed at default settings.

Data were presented as mean ± SEM of independent biological replicates, unless otherwise noted. Distribution patterns of the samples were evaluated through Shapiro-Wilk test (Shapiro and Wilk, 1965) and those populations that followed normal distribution were analyzed by two-tailed unpaired Student’s t-test. Samples that did not pass normality test were analyzed by nonparametric two-tailed Mann-Whitney U test (Mann and Whitney, 1947). p-values of ≤0.05 were considered as significantly different. All analyses were performed through GraphPad Prism 8 software and all values are displayed in Supplementary Data S3.

In order to query for Ae. albopictus lncRNAs that are potentially implicated in mosquito reproduction, we initiated our computational analysis on the respective NCBI gene models, which contain 8571 transcripts that are predicted to belong to the ncRNA class (“gbkey = ncRNA”). Bioinformatics approaches were adopted to annotate further these RNAs regarding their coding potential, species specificity, genomic localization, as well as their context-expression specificity within Ae. albopictus.

Coding potential estimates of the predicted ncRNAs were calculated with CPC2 and FEELnc. CPC2 (Kang et al., 2017) is a well-known species-agnostic Support Vector Machine model that makes use of Fickett score, open reading frame (ORF) length and integrity and isoelectric point to predict coding potential. In contrast, FEELnc (Wucher et al., 2017) accepts as input known coding and non-coding transcripts from a species of interest and trains a species-aware coding potential model (Random Forests) using as predictors the ORF coverage of transcripts, k-mer composition (1-, 2-, 3-, 6-, 9-, and 12-mers) and transcript length. When known non-coding RNA sequences are not available, under shuffle mode, FEELnc shuffles the provided protein coding sequences but preserves 7-mer frequencies, to create a faux set of non-coding transcripts to use in model training-testing, a method shown to fare adequately in benchmarks against real lncRNAs (Wucher et al., 2017). The FEELnc model output is a score between 0 (i.e., no coding potential) and 1 (coding potential), while its optimal cutoff point that maximizes sensitivity and specificity is calculated by means of 10-fold cross validation. CPC2 annotated 8289 “ncRNA” transcripts as non-coding and 282 as potentially coding. FEELnc classification yielded 8096 and 475 transcripts labeled as non-coding and coding respectively, at an optimal cutoff point of 0.43 (Figure 1A). CPC2 and FEELnc predictions were juxtaposed with each other, highlighting 7,900 and 86 transcripts (92 and 1% of total) concordantly classified as non-coding and coding RNA respectively, and identifying 589 predictions (7% of total) that were discrepant between the two approaches (Figure 1B, Supplementary Data S1).

In order to identify species-specific lncRNAs, transcript homology queries were performed (local BLASTn) against all Hexapoda Reference RNA sequences. Each Ae. albopictus “ncRNA” transcript was annotated with respect to the number of other species it presented hits against. Out of the total transcripts, 7,827 (91%) were found to exhibit no similarity with RNAs of any other Hexapoda species (Figure 1C).

FEELnc was also utilized to annotate transcripts’ localization with regard to transcripts annotated as protein coding (Figure 1D). In total, 2051 instances (24%) were found to overlap protein coding transcripts in sense (9%) or antisense (15%) orientation. The remaining 76% was divided among intergenic transcripts that presented non-overlapping neighboring genes within 100 kb of their loci (5,318 transcripts, 62%) and those that were found to exist in genomic regions not harboring other genes (annotated as “distal” intergenic, 1,202 instances, 14%).

Finally, in order to obtain expression metrics with regard to these ncRNAs in discrete developmental stages and tissues, a publicly available developmental transcriptomic dataset of Ae. albopictus (Gamez et al., 2020) was analyzed from scratch (SRA accession number SRP219966). Briefly, the dataset captured expression estimates from adult ovaries of non-blood-fed (NBF, fed with a 10%-sucrose solution) and post-blood meal (PBM, at 12-24-36-48-60-72 time-points) insects, carcasses (i.e., female bodies without the ovaries, also NBF and PBM at the same time-points), adult testes, diapause, larvae, pupae and embryo samples at numerous time-points. In the absence of replicates (with the exception of testes which were in duplicate), no attempt to assess differential expression was performed. Instead, log₂-transformed TMM-normalized counts were used to obtain tau indices (Yanai et al., 2005), which have been recently shown to perform consistently as a tissue specificity metric (Kryuchkova-Mostacci and Robinson-Rechavi, 2017). Briefly, tau constitutes a per-transcript tissue specificity unit ranging from 0 (no specificity) to 1 (highest specificity), while the fractions of each tissue that contributed to calculation of tau can also be obtained. Within the current dataset that is composed of a number of distinct tissues, developmental stages and time-points, we denoted that tau be regarded as a context-specificity index. The fractions of 3,119 transcripts exhibiting tau > 0.5 are presented as a heatmap (Figure 1E) for all available contexts, outlining the existence of numerous instances of intermediate-to-high specificity in ovary contexts (i.e., yellow-to-red transcripts within ovary samples).

We focused our attention on 984 transcripts that were 1) annotated as having no/limited coding potential by both CPC2 and FEELnc, 2) presented no sequence similarity to known transcripts of other Hexapoda and 3) exhibited tau > 0.5 and tau fraction > 0.5 in at least one PBM ovary time-point. Aiming to experimentally scrutinize a limited number of lncRNAs, out of this subset we selected 10 ovary-specific lncRNAs that were upregulated upon blood-feeding and presented limited or no expression in other developmental samples (Supplementary Table S1; Figure 2A). We termed these lncRNAs Norma (NOn-coding RNAs in Mosquito ovAries).

In order to clarify the potential role of the 10 Norma genes, we evaluated their phenotypic impact in reproduction through RNAi silencing. Given their expression pattern, we reasoned that knocking down Norma genes would mostly impact mosquito oviposition. Τo assess this, we generated in vitro dsRNA against each of the ten Norma lncRNAs which we injected into 5-day old inseminated adult females. We then provided a blood meal to the injected mosquitoes and monitored them for oviposition. We counted the number of eggs that were laid by each individual mosquito treated with any of the ten Norma-dsRNA (anti-Norma) and compared oviposition with a GFP-dsRNA sample (anti-GFP). A non-injected, blood-fed sample (untreated) was also present in the assay to monitor the effect of the environmental conditions. We observed that females injected with six different Norma-dsRNAs laid fewer eggs compared to GFP-dsRNA injected control, indicating the potential influence of silencing the corresponding lncRNAs to reproduction (Figure 2B). Specifically, dsRNA against Norma3, Norma9, Norma16, Norma17, Norma18 & Norma19 exhibited statistically significant reduction of oviposition rates compared to the GFP control (p-value ≤0.05). We focused our downstream analyses on Norma3 because it also presented some further striking features: an ovary-specific expression profile along with a sharp upregulation pattern in the post-vitellogenic time-points. Other lncRNAs that influenced oviposition (Norma9, Norma16, Norma17, Norma18, and Norma19) will be the subject of a future investigation.

Initially, we determined the detailed expression profile of Norma3 by examining its stage- and tissue-related patterns in samples collected from our laboratory mosquito strain. We collected ovaries at the same time-points as the ones that were described in the RNA seq dataset that we analyzed. Specifically, we collected ovaries from NBF and PBM (12, 24, 36, 48, 60, and 72 h) time-points. In addition, we dissected and stored other tissues from the same time-points (carcass, midgut, Malpighian tubes, and head). Norma3 exhibited basal expression in NBF and early PBM samples, while its expression abruptly increased >2000-fold in 60 h-PBM ovaries, compared to the expression in NBF ovaries, and then significantly dropped to 200-fold in 72 h-PBM ovaries, compared to NBF (Figure 3A). At the same time-point (60 h PBM), Norma3 displayed basal expression in all other examined tissues (carcass, midgut, Malpighian tubes, head) (Figure 3A). The expression pattern was in accordance with publicly available RNA seq data (Gamez et al., 2020). Subsequently, we estimated the efficiency of Norma3 silencing in the ovaries of female mosquitoes injected with Norma3-dsRNA (anti-Norma3), compared with mosquitoes injected with GFP-dsRNA (anti-GFP). We collected ovaries from individual mosquitoes injected with either Norma3-dsRNA or GFP-dsRNA 60 h PBM and measured the effect of RNAi on the expression levels of Norma3. Results showed an average expression drop of >50% in anti-Norma3 replicates compared to anti-GFP (Figure 3B), which presented adequate statistical significance (Student’s t-test, p = 0.0031) to support our experimental pipeline.

To more deeply characterize the impact of Norma3 silencing on Ae. albopictus reproductive ability, we examined various phenotypic traits that are connected to reproduction. First, we looked at the ovary morphology at 60 h PBM by measuring the long axis of the ovoid follicle of ovaries dissected by individuals of the anti-Norma3 dsRNA and anti-GFP dsRNA samples. We repeated the dissection process three times and detected that ovaries of anti-GFP sample presented an average length of 293.6 nm, while ovaries of anti-Norma3 presented a significantly lower length of 240 nm (Mann-Whitney test, p = 0.003) (Figure 4A). Interestingly, anti-Norma3 sample included several much smaller ovaries that had a shorter follicular diameter (Supplementary Data S3), while there was evidence of the presence of nurse cells along with their oocyte, indicating their delayed development (Figure 4A,B, Supplementary Figure S1B). Nine anti-Norma3 smaller ovaries that displayed a mean length of 153.4 nm were collected and stored for further analysis.

Afterwards, we estimated the effect of the Norma3-dsRNA on oviposition by counting the number of eggs that were laid individually by each female of anti-Norma3, anti-GFP and untreated samples. Mosquitoes of the untreated control laid an average of 44.6 eggs, while mosquitoes of the anti-GFP control laid an average of 41.6 eggs. On the other hand, mosquitoes of the anti-Norma3 treatment laid an average of 24.1 eggs. Interestingly, 14% (n = 15) of the anti-Norma3 replicates laid zero eggs, while 26% (n = 27) laid fewer than 10 eggs. In contrast, 3% (n = 3) of the anti-GFP and 1% of the untreated samples (n = 1) laid 0 eggs, while 6% (n = 6) of the anti-GFP and 3% of the untreated samples (n = 3) laid fewer than 10 eggs (Supplementary Data S3). The reduction of the average oviposition rate between anti-Norma3 and anti-GFP was 43% and exhibited high statistical significance (Student’s t-test, p < 0.0001) (Figure 4C).

Subsequently, we addressed the hatch rate of untreated, anti-GFP and anti-Norma3 samples. We counted the number of larvae that emerged from the eggs that were laid by each individual mosquito of the studied samples and we divided by the total amount of eggs that were laid by each mosquito. Untreated mosquitoes displayed an average hatch rate of 79%, anti-GFP mosquitoes presented a similar hatch rate 78%, while anti-Norma3 mosquitoes had a much lower rate of 37%. (Figure 5A). The difference of the hatch rate between anti-GFP and anti-Norma3 was about 53% and presented high statistical significance (Mann-Whitney test, p < 0.0001).

Finally, we attempted to visualize possible larval defects in the anti-Norma3 treated sample. For this, we rendered 72 h eggs transparent by bleaching and observed the larvae under the microscope. We investigated the entire batches of eggs/embryos laid by each mosquito treated with dsRNA either against Norma3 or GFP and detected significant changes. There was a constant percentage of 75-85% that presented the deficient phenotype in anti-Norma3 sample, while the respective percentage of anti-GFP embryos was 10-15%. We dechorionated the embryos laid by 10 anti-Norma3 and 10 anti-GFP mosquitoes to validate the result. Representative images of dechorionated embryos that were laid by different mosquitoes of each sample and displayed the characteristic phenotype are presented in Figure 5B and Supplementary Figure S2. In anti-GFP embryos the eight abdominal segments, the thoracic segments, the respiratory siphon, the head and the ocelli were clearly visible indicating regular development of the embryo. On the other hand, anti-Norma3 embryos did not present any of the anticipated normal structures; instead, they exhibited a defective appearance of an undifferentiated mass (Figure 5B, Supplementary Figure S2B).

Since lncRNAs often regulate coding genes in their genomic neighborhood (Rinn and Chang, 2012; Engreitz et al., 2016; Joung et al., 2017; Khyzha et al., 2019; Xing et al., 2021), we investigated the possible association of Norma3 with protein-coding genes in its vicinity. We returned to the available transcriptomic study containing NBF and PBM ovary datasets (Gamez et al., 2020) and assessed the correlation of Norma3 expression against the rest of the transcriptome. Within the region from 100 kb upstream to 100 kb downstream of Norma3, we identified 8 annotated transcripts (i.e., 4 mucins, 1 venom protein, and 3 chymotrypsin inhibitors) which all exhibited robust positive correlation with Norma3 (Pearson correlation coefficient >0.96, maximum FDR = 2.17e-9) (Supplementary Data S2). Subsequently, we reduced the transcriptome to transcripts which were more likely to present any change in ovaries among early (i.e., 12–24 h), middle (36–48 h) and late (60–72 h) time points (posterior-probability EBSeq >75%). This transcript set (n = 10,314) was subjected to clustering analysis based on their expression over the entire time course (i.e., NBF, 12, 24, 36, 48, 60, and 72 h PBM). Norma3 was grouped together with 114 other genes (111 protein-coding and 3 long non-coding) in a cluster of genes that was not expressed in NBF and early PBM samples but initiated low transcription at 36 h and peaked at the end of vitellogenesis (around 48 h PBM) (Figure 6A). Among the 8 positively correlated neighbor transcripts, three mucins belonged to a cluster that presented an intense expression peak at 48 h (Figure 6B).

We then generated their detailed expression profile in NBF and PBM ovaries. We confirmed the expression pattern of five of the neighboring transcripts via qPCR. Three were annotated as mucin-2 and mucin-2-like (mucin1, 2, 3), one was annotated as cysteine-rich venom protein 6-like (venom1) and one as chymotrypsin inhibitor (chymo1). All five genes presented basal expression in NBF, 12, 24, 36, and 48 h ovaries. The expression of the three mucins presented a sharp increase between 12,000 (mucin2) and ∼40,000-fold (mucin3) in 60 h PBM ovaries, while chymo1 increased ∼10,000-fold in 72 h PBM and venom1 3,000-fold in 60 h PBM ovaries (Supplementary Figure S3). To further assess the potential impact of Norma3 on the mucins we assayed the expression of the three mucins in a small sample of smaller ovaries (Figure 4A) of the anti-Norma3 60 h-PBM sample (that we collected as described in Section 3.4). Statistically significant downregulation between 40 to 80% of two out of the three mucin genes was detected in the smaller anti-Norma3 ovaries compared to anti-GFP control (Figure 6C). Downregulation of mucin3 presented the most significant effect, while no significant effect was detected on other protein-coding transcripts located in the same genomic region (mucin1, chymo1, venom1). Nevertheless, since we only tested nine smaller anti-Norma3 ovaries, we acknowledge that this result serves as a preliminary indication of the potential impact of Norma3 on mucin expression that should be validated with larger datasets.