Post-transcriptional regulation of frog innate immunity: discovery of frog microRNAs associated with antiviral responses and ranavirus infection using a Xenopus laevis skin epithelial-like cell line

Xela DS2, a skin epithelial-like cell line previously generated from X. laevis dorsal skin (Bui-Marinos et al. 2020), was maintained in seven parts Leibovitz’s L-15 medium (Wisent, Mont-Saint-Hilaire, Canada) diluted with 3 parts sterile cell culture water to adjust for amphibian osmolarity, and is herein referred to as amphibian L-15 (AL-15) medium. AL-15 medium was supplemented with 15% fetal bovine serum (FBS; lot #234K18; VWR, Radnor, United States) and used to culture Xela DS2 cells at 26 °C in 75 cm² plug-seal tissue-culture treated flasks (BioLite, Thermo Fisher Scientific, Waltham, United States). Cells were split 1:4 upon reaching ∼90%–95% confluency, which usually occurred every 3–4 d. Epithelioma Papulosum Cyprini (EPC) cells were maintained in Leibovitz’s L-15 medium supplemented with 10% FBS. EPC cells were cultured in 75 cm² plug-seal tissue-culture treated flasks at 26 °C and were split 1:4 weekly. Prior to use in experiments, viable Xela DS2 and EPC cells were enumerated using Trypan blue (0.2% final concentration; Invitrogen, Waltham, United States) and a haemocytometer, and seeded to account for their respective plating efficiencies (79% for Xela DS2 and 100% for EPC).

FV3 (Granoff strain) was propagated in EPC cells as described previously (Bui-Marinos et al. 2021). Briefly, a monolayer of EPC cells was infected with a 1:10 dilution of stock FV3 in L-15 medium supplemented with 2% FBS for seven days at 26 °C. Seven days post-infection, virus-containing cell culture medium was collected, centrifuged at 1000 × g for 10 min, and filtered through a 0.22 μm PES filter (FroggaBio, Concord, Canada) prior to storage at −80 °C until use.

The tissue culture infectious dose values, wherein 50% of cells are infected (TCID₅₀/mL), were determined according to the modified Kärber method (Kärber 1931; Pham et al. 2011). EPC cells (100 000/well) were seeded in a 96-well plate (BioBasic, Toronto, Canada) and allowed to adhere overnight at 26 °C. The following day, cell culture medium was removed, and cells were treated with 200 μL of a 10-fold dilution series of the virus-containing cell culture supernatants to be tested (FV3 stock prepared as described above or experimental sample supernatant) diluted in L-15 medium containing 2% FBS. Plates were sealed with parafilm and incubated at 26 °C for seven days prior to visual examination of cytopathic effects (CPE). EPC monoloayers were scored for CPE and approximate FV3 plaque forming unit (PFU)/mL values were calculated by multiplying TCID₅₀/mL values by a factor of 0.7 (Knudson and Tinsley 1974). PFU/mL values were used to calculate multiplicity of infection (MOI) (Knudson and Tinsley 1974). The FV3 stock used in this study was determined to have a TCID₅₀/mL value of 1.64 × 10⁹.

Xela DS2 cells (passages 77, 80, and 83; n = 3) were seeded into two 6-well plates (700 000 cells/well) in AL-15 medium supplemented with 15% FBS and allowed to form a monolayer overnight at 26 °C. The following day, cell culture medium was removed, and cells were treated with 2 mL of medium alone (AL-15 supplemented with 2% FBS, two wells per plate) or medium containing FV3 at an MOI of 2 (FV3 infected, one well per plate). The amount of virus required to obtain an MOI of 2 was calculated by first converting TCID₅₀/mL values to PFU/mL values (Knudson and Tinsley 1974). The raw TCID₅₀ values added to the Xela DS2 cells, without conversion to PFUs, was 2.0 × 10⁶. After 2 h of incubation at 26 °C, the medium was removed, cells were washed three times with 2 mL amphibian phosphate buffered saline (8 g/L sodium chloride, 0.5 g/L potassium chloride, 2.68 g/L sodium phosphate dibasic heptahydrate, 0.24 g/L potassium phosphate monobasic), and 2 mL of fresh AL-15 medium supplemented with 2% FBS was added to all wells. At this point, 1 μg/mL of poly(I:C) composed of low and high molecular weight molecules (cat # P1530; Sigma-Aldrich, Burlington, United States) was added to one of the untreated wells. Untreated, poly(I:C)-treated, and FV3-infected cells were incubated at 26 °C for 24 h, while another set of untreated and FV3-infected cells were incubated at 26 °C for 72 h. Prior to collection of medium to determine TCID₅₀/mL values and cells for RNA extraction, cell images were captured using a Leica DMi1 phase-contrast microscope fitted with an MC170 colour camera and LASX 4.8 software.

Adherent Xela DS2 cells and any detached cells in suspension were harvested for total RNA isolation using the Monarch Total RNA Minipreps kit (New England Biolabs, Ipswich, United States) according to the manufacturer’s recommended protocol. Suspension cells present in the cell culture medium were collected by centrifugation at 500 × g for 2 min and the resulting cell culture medium was aspirated. Lysis buffer (800 μL) was added to each well and incubated with adherent cells for 2 min at room temperature. Lysed cells were transferred from the 6-well plates to microcentrifuge tubes containing the corresponding pelleted suspension cells. Cells were resuspended by pipetting and incubated for 2 min to permit cell lysis. At this stage, an exogenous miRNA (cel-miR-39 Spike-In kit; Norgen Biotek Corp., Thorold, Canada) was spiked in (99 fmol) to serve as a positive control for miRNA detection. On-column DNase I treatment was performed for enzymatic removal of residual genomic DNA according to the manufacturer’s instructions, and elution was performed in nuclease-free water (50 μL). Total RNA was quantified using the Take3 microvolume plate accessory on a Cytation 5 Cell Imaging Multi-Mode Reader (Biotek, Winooski, United States), RNA purity was determined by analyzing A_260/280 and A_260/230 ratios, and RNA integrity was examined by electrophoresis on a bleach agarose gel (Aranda et al. 2012).

For each sample, ∼2 μg of RNA (∼80 ng/μL in 25 μL) was shipped to The Centre for Applied Genomics (TCAG) at The Hospital for Sick Children (Toronto, Canada) on dry ice. TCAG performed further quality assessment using a Bioanalyzer 2100 and RNA 6000 Nano LabChip Kit (Agilent Technologies Inc., Santa Clara, United States), and samples with an RNA integrity number >8.5 were used for small RNA-seq library preparation. Library preparation and sequencing was performed by TCAG. Small RNA-seq libraries were prepared using the NEBNext Multiplex Small RNA Library Prep Set for Illumina (New England Biolabs, Ipswich, United States). Small RNAs were enriched using bead-based size selection. Resulting small RNA libraries were sequenced on an Illumina HiSeq 2500 rapid-run flow cell to generate 50 nucleotide (nt) single-end reads.

Raw ∼50 nt single-end sequencing reads were adaptor-trimmed and size-filtered using Cutadapt (Martin 2011) with the following parameters: single-end reads, 3′ adaptor of AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC, minimum length of 17, maximum length of 30. Read quality reports were generated using FastQC (Andrews 2010) with default parameters. High quality (Q > 34; 0.04% error rate) clean reads were mapped to the X. laevis genome (J strain; GCF_001663975.1) using Bowtie2 (Langmead and Salzberg 2012) with the following parameters: single-end reads, write aligned and unaligned reads to separate files, allow one mismatch, output in sequence alignment map (SAM) format. Reads that mapped to the X. laevis genome were retained for further analyses. Cutadapt, FastQC, and Bowtie2 were run using the Galaxy web platform (usegalaxy.org; (Afgan et al. 2016)).

Novel X. laevis miRNAs were detected using miRDeep2 v0.1.3 (Friedländer et al. 2012). Aligned reads were collapsed into unique sequences using the mapper.pl script and reads with counts <10 in a given library were discarded for this purpose. Transcripts with read counts <10 are generally considered to be noise and are often removed prior to downstream bioinformatics analysis (Law et al. 2016); thus, removing these reads prior to novel miRNA prediction helps improve the signal-to-noise ratio, which is an important component of miRDeep2-based miRNA detection. The miRDeep2.pl script was executed using default parameters and the following inputs: collapsed reads, alignments, the X. laevis genome, known X. laevis pre-miRNAs (miRBase), and known X. laevis and Xenopus tropicalis (Western clawed frog) miRNAs (miRBase). Secondary structures associated with novel X. laevis pre-miRNAs were determined using the RNAFold component of the ViennaRNA package v2.4.17 (Lorenz et al. 2011) under default parameters. Novel X. laevis miRNAs were considered high confidence if they were characterized by a miRDeep2 score ≥ 4 (≥ 90% probability of being a true positive) or if their sequence perfectly matched a known X. tropicalis miRNA.

X. laevis miRNA read counts were generated using miRDeep2 (Friedländer et al. 2012). Collapsed X. laevis reads, known and novel X. laevis miRNA sequences, and X. laevis pre-miRNA sequences (known (retrieved from miRBase) and novel (predicted by miRDeep2)) were input into the miRDeep2 quantifier.pl script using default parameters. Reads corresponding to several miRNAs were in very low relative abundance (e.g., <10 in a given library). To improve the accuracy of our analyses by strengthening the signal-to-noise ratio (Law et al. 2016), miRNAs were considered detected if their read counts were at least 10 in at least one library. miRNAs with <10 read counts in all libraries were considered undetected and were excluded from downstream analyses. Read counts served as input for differential expression analysis using DESeq2 (Love et al. 2014) and EdgeR (Robinson et al. 2010). DESeq2 parameters included “select datasets per level”, input data = count data, and visualizing results = yes. EdgeR parameters included “single count matrix”, use factor information file = yes, use gene annotation = no, normalization method = TMM. Known and novel X. laevis miRNAs were considered differentially expressed if a false discovery rate (FDR) <0.05 was reached by both programs (consensus DE). Heatmaps were generated using heatmap2 with the following parameters: enable data clustering = no, coloring groups = blue to white to red. DESeq2, EdgeR, and heatmap2 were run using the Galaxy web platform (usegalaxy.org; (Afgan et al. 2016)).

X. laevis 3′ UTRs were retrieved from the University of California Santa Cruz (UCSC) Table Browser using the Xenbase trackhub (genome.ucsc.edu; genome assembly v9.2). FV3 coding sequences (cds) were retrieved from National Center for Biotechnology Information (NCBI) (Accession NC_005946.1). Prior to running target prediction analysis, the dinucleotide frequencies of X. laevis 3′ UTRs and FV3 cds were calculated using the compseq function of EMBOSS v6.6.0 (Rice et al. 2000). Using these frequencies in target prediction analyses improves p-value accuracy (Rehmsmeier et al. 2004). The targets of known and novel X. laevis miRNAs were predicted using RNAcalibrate and RNAhybrid v2.1.2 (Rehmsmeier et al. 2004). RNAcalibrate was first used to estimate distribution parameters using default parameters (in addition to strict seed pairing at bases 2–7 of the miRNA and the dinucleotide frequencies generated by EMBOSS). Next, RNAhybrid was used to predict the X. laevis and FV3 mRNA targets of known and novel X. laevis miRNAs. RNAhybrid parameters included maximum energy = −25 kcal/mol, maximum mismatches in internal loop = 1, strict seed pairing at miRNA bases 2–7, mismatches in bulge loop = 0, one hit per interaction, p < 0.01. Gene functions were annotated using Xenbase (Karimi et al. 2018). Immune-related X. laevis targets were identified through a combination of manual search and cross-referencing The Gene Ontology (GO) Resource (The Gene Ontology Consortium 2019) and InnateDB (Breuer et al. 2013).

The PANTHER overrepresentation test was used to perform GO term enrichment analysis to reveal biological processes that are overrepresented in the list of X. laevis genes predicted to be targeted by differentially expressed X. laevis miRNAs (The Gene Ontology Consortium 2019). PANTHER parameters were as follows: reference list: all genes in database, “GO biological process complete” annotation dataset, Fisher’s exact test, FDR <0.05. GO term analysis was performed at geneontology.org.

RT-qPCR experiments were performed on a duplicate set of RNA samples isolated from cells treated in parallel with RNA-seq samples. Total RNA (1 μg) was reverse-transcribed using the microScript miRNA cDNA synthesis kit (Norgen Biotech Corp., Thorold, Canada) according to the manufacturer’s recommended protocol. Resulting cDNA was diluted (1/20) with nuclease-free water and stored at −20 °C until RT-qPCR analysis. RT-qPCR reactions consisted of 1 × PowerUP SYBR green mix (Thermo Fisher Scientific, Waltham, United States), 2.5 μL of diluted cDNA, and 500 nM (each) forward (miRNA-specific) and reverse (universal) primers in a final reaction volume of 10 μL. RT-qPCR primer sequences and efficiencies are presented in Supplementary Table S1. All RT-qPCR reactions were performed in triplicate on a QuantStudio5 Real-Time PCR System (Thermo Fisher Scientific, Waltham, United States), and reaction conditions were as follows: denaturation at 50 °C for 2 min and 95 °C for 2 min, followed by 40 cycles of denaturation at 95 °C for 1 s and extension at 60 °C for 30 s. A melt curve step (denaturation at 95 °C for 1 s and dissociation analysis at 60 °C for 20 s followed by 0.1 °C increments from 60 °C to 95 °C at 0.1 °C/s) was used to ensure the presence of a single dissociation peak.

RT-qPCR data were analyzed using the ΔΔCt method in Microsoft Excel with xla-miR-16a-5p serving as an endogenous control. xla-miR-16a-5p was chosen as a control for RT-qPCR normalization by examination of the small RNA-seq data (normalized reads per million (RPM) read counts generated by the miRDeep2 quantifier.pl script). This control miRNA was chosen based on the following criteria: (i) <10% variability in expression (relative to total number of host-mapped reads) between experimental libraries and their respective control libraries and (ii) >1000 RPM in the libraries of at least one treatment. xla-miR-16a-5p served as the endogenous control for RT-qPCRs targeting miRNAs, as well as dicer1 and drosha transcripts. While we recognize that normalization of an mRNA to a miRNA control is not ideal, it is difficult to identify an appropriate endogenous control gene for use in infection experiments, particularly late in viral infection. As the relative proportion of host RNA in total RNA samples decreases as viral RNA levels increase in infected cells, it is difficult to experimentally validate an endogenous control gene by RT-qPCR alone, due to the inherent instability of host RNA abundance between untreated and infected cells. Therefore, small RNA-seq data was used to determine an appropriate endogenous control miRNA in lieu of an mRNA control gene, which permitted accurate assessment of the stability of the control miRNA relative to the total host RNA abundance in each sample. Statistically significant differences were identified using one-way ANOVA tests paired with Dunnett’s post-hoc tests or Student’s t-tests where appropriate, and differences were considered statistically significant if p < 0.05.

The small RNA-seq datasets generated by this study are available in the NCBI Short Read Archive repository (BioProject ID #PRJNA752644).

To identify miRNAs that are involved in antiviral responses in frog skin epithelial cells, small RNAs from untreated (24 h and 72 h), poly(I:C)-treated (24 h), and FV3-infected (24 h and 72 h) Xela DS2 skin epithelial-like cells were sequenced. A poly(I:C) concentration of 1 μg/mL and a treatment duration of 24 h were chosen as these conditions produce robust transcriptional immune responses in Xela DS2 cells (Bui-Marinos et al. 2020), pretreatment with poly(I:C) confers functional protection against FV3 infection in Xela DS2 (Bui-Marinos et al. 2021), and similar concentrations have been used in previous studies profiling miRNA responses to poly(I:C) treatment in other species (Wang et al. 2017; Wu et al. 2019). Poly(I:C), comprised of a mixture of low molecular weight and high molecular weight molecules, was used to roughly mimic the viral RNAs of various lengths produced during FV3 replication (Tan et al. 2004). An FV3 MOI of 2 was chosen to approach a synchronous infection of Xela DS2 (approximately 86% of cells should receive one or more viral particles as calculated using the Poisson distribution), and thus facilitate the detection of host miRNAs in response to viral infection while minimizing potential cytotoxic effects of higher doses of the virus. In the present study, we opted not to examine FV3 infection of Xela DS2 cells at a lower MOI as that would result in an asynchronous infection and may obscure information that could be obtained regarding the timing of cellular miRNA expression in relation to progression of viral replication within host cells. Sampling time points of 24 h and 72 h were chosen to capture differences in antiviral responses at different times post infection as overt FV3-induced CPE are not obvious at earlier time points (e.g., 12 h), moderate CPE is observed at 24 h, and extensive CPE is evident at 72 h in Xela DS2 cells infected with FV3 at an MOI of 2 (Supplementary Fig. S1a). A later time point (i.e., 96 h) was not included based on previous observations of loss of cell viability in FV3-infected cells (Bui-Marinos et al. 2021). The TCID₅₀/mL values measured following FV3 infection of Xela DS2 cells at an MOI of 2 for 24 h or 72 h were 6.95 × 10⁵ (log₁₀ = 5.84) or 2.81 × 10⁷ (log₁₀ = 7.44), respectively (Supplementary Fig. S1b). Small RNAs were sequenced on an Illumina HiSeq 2500 resulting in 15 small RNA libraries (n = 3 per treatment; Supplementary Fig. S1c–d). A schematic of the bioinformatics pipeline employed in this study is depicted in Fig. 1a. Approximately 7.6–11.3 million raw reads per sample underwent adapter, length, and quality filtering using Cutadapt (Martin 2011) and FastQC (Andrews 2010) to obtain ∼1.5–5.9 million clean reads (17–30 nt) per library (Fig. 1b). Genome alignment was performed using Bowtie2 (Langmead and Salzberg 2012). Reads that aligned to the X. laevis genome were subjected to further analyses. The number of X. laevis-mapped reads in each library ranged from ∼29% to 97% (Fig. 1b), and their size distribution spanned from 17 nt to 30 nt (Fig. 1c). More than 79% of the small RNA reads that mapped to the X. laevis genome were 21–24 nt in length, with 22–23 nt RNAs representing the majority.

At present, there are 247 mature X. laevis miRNAs in the miRBase database (Kozomara et al. 2019), a number which is relatively low for vertebrates. Thus, we sought to first identify novel X. laevis miRNAs from the small RNA libraries. Using miRDeep2 software, 133 novel X. laevis miRNAs were detected, which were classified into 52 families (Supplementary Table S2). Twenty-five of these novel X. laevis miRNAs are conserved in other species and present in miRBase, while the remainder are entirely novel miRNAs in any species (Fig. 1d). Many of these entirely novel miRNAs are derived from repetitive genomic regions (Supplementary Table S2). Conserved miRNAs were assigned the name of the corresponding orthologous miRNA in miRBase (Kozomara et al. 2019), while truly novel miRNAs were numbered sequentially and given the prefix “novel”. There is a bias towards “U” as the first (5′) nucleotide of these novel X. laevis miRNAs, which is characteristic of vertebrate miRNAs (Yang et al. 2011). Next, miRDeep2 software was used to detect and quantify known and novel X. laevis miRNAs in each library. Along with the 133 novel X. laevis miRNAs, 136 known X. laevis miRNAs were detected across all 15 libraries (Fig. 1d). Thus, a total of 269 X. laevis miRNAs were detected in the small RNA-seq datasets. Most of these miRNAs were shared amongst all groups; however, two miRNAs were uniquely detected in FV3-infected cells at 24 h and 72 h, while one miRNA was uniquely detected in untreated cells at 72 h (Fig. 1e).

Pathogens are sensed by a variety of PRRs that are either cytosolic (e.g., cGAS, RIG-I/MDA5) or membrane-associated (e.g., toll-like receptors (TLRs)) and recognize pathogen-associated molecular patterns (PAMPs) such as viral nucleic acids (Mogensen 2009; Motwani et al. 2019). PAMP sensing by PRRs can initiate antiviral signaling cascades that culminate in the production of IFNs and inflammatory cytokines that combat infection (Mogensen 2009). To gain a clearer picture of the antiviral pathways targeted by miRNAs that are differentially expressed in response to poly(I:C) or FV3, we mapped these miRNAs onto several key antiviral pathways such as cGAS-STING, RIG-I/MDA5, TLR3, TLR4, TLR7, TLR8, TLR9, and type I IFN signaling pathways (Fig. 4). Poly(I:C)-altered miRNAs target components of all eight of these pathways (stim1, trim25, tbk1, tram1, ripk1, tab2, irak2, irf5, nmi, ifnar1, stat1, ubr4, and pkc-δ), as well as their downstream products such as NF-ĸB-induced genes and ISGs (Fig. 4a). On the other hand, FV3-altered miRNAs targeted fewer components of these pathways (trim14, stim1, trim25, tbk1, tram1, ripk1, tab2, stat1, ubr4, and pkc-δ) and fewer NF-ĸB-induced genes and ISGs (Fig. 4b). While some targets in these pathways are found in common between miRNAs altered by poly(I:C) and FV3, there are several distinct differences. Poly(I:C)-altered miRNAs uniquely target irak2, irf5, nmi, and ifnar1, while FV3-altered miRNAs uniquely target trim14. Interestingly, while several poly(I:C)-altered miRNAs target components of TLR7/8/9 signaling, FV3 is associated with a lack of changes in the expression of miRNAs targeting these pathways.

RT-qPCR was used to validate the differential expression profiles of a subset of X. laevis miRNAs that are predicted to target important regulators of innate antiviral responses. Twelve miRNAs were chosen that (i) are predicted to target important immune-related genes; (ii) represent a combination of known, conserved, and novel miRNAs; (iii) represent a combination of upregulated and downregulated miRNAs; and (iv) represent a combination of miRNAs with changes in expression that are consistent between the poly(I:C) and FV3 treatment groups or are inconsistent between these groups. These miRNAs included xla-miR-106-5p, xla-miR-130c-5p, xla-miR-139-3p, xla-miR-140-5p, xla-miR-222-3p, novel-xla-miR-7-5p, xla-miR-22-5p, novel-xla-miR-12-3p, novel-xla-miR-16-5p, novel-xla-miR-16-3p, novel-xla-miR-28-5p, and novel-xla-miR-28-3p. The RNA-seq-derived differential expression profiles of 11 of these miRNAs were successfully validated by RT-qPCR (Fig. 5), and one miRNA (xla-miR-222-3p) did not show comparable differential expression profiles between the two techniques. With regards to the 11 validated profiles, there are a few cases where one technique produced results that were statistically significant whereas the other technique did not (xla-miR-106-5p in FV3-infected cells at 24 h and 72 h; novel-xla-miR-12-3p in FV3-infected cells at 72 h, novel-xla-miR-16-5p in FV3-infected cells at 24 h, and novel-xla-miR-16-3p in FV3-infected cells at 24 h). These are likely due to slight variation in replicates and differences in sensitivity of the two techniques. For example, novel-xla-miR-16-5p and novel-xla-miR-16-3p were undetected in untreated samples by RNA-seq, whereas amplification was achieved using RT-qPCR. The lack of detection in RNA-seq can slightly skew the fold-change values calculated. Despite the sometimes slight variance in magnitude of differential expression between the two techniques (RT-qPCR and RNA-seq), the overall trends remained the same.

miRNAs that target immune-related genes were more frequently downregulated than upregulated in response to either poly(I:C) or FV3. In line with this, and as observed in other pathosystems, host cells have been found to modulate global miRNA biogenesis to manipulate miRNA expression during infection (Bruscella et al. 2017). Thus, we sought to profile changes in the transcript levels of components of the miRNA biogenesis pathway, dicer1 and drosha, in response to poly(I:C) or FV3. RT-qPCR analysis revealed downregulation of dicer1 transcripts in response to poly(I:C) (24 h) as well as FV3 (24 and 72 h) (Fig. 6a). In addition, drosha transcripts were downregulated in response to FV3 infection at 72 h (Fig. 6b). Together, these results suggest that miRNA biogenesis is partially suppressed in response to viral infection.

We uncovered 80 X. laevis miRNAs that are differentially expressed in response to extracellular poly(I:C) and (or) FV3 in a frog skin epithelial-like cell line, suggesting for the first time that frog miRNAs are key regulators of antiviral immunity. Changes in miRNA expression in response to poly(I:C) were profiled to model robust antiviral responses, as poly(I:C) induces potent antiviral responses and is an inert molecule (nonimmunoevasive). Extracellular delivery of poly(I:C) is intended to mimic the viral dsRNA replicative products released into the extracellular space following lysis of infected cells and is subsequently sensed by neighboring cells. Poly(I:C) treatment of Xela DS2 cells altered the levels of 65 X. laevis miRNAs, which may play a role in post-transcriptional regulation of host gene expression in response to exogenous viral dsRNA. Poly(I:C) treatment of Xela DS2 induced changes in miRNAs that are predicted to target important regulators of antiviral signaling pathways such as cGAS-STING (e.g., stim1), RIG-I/MDA5 (e.g., trim25), TLR signaling (e.g., irak2, irf5, ripk1, and tab2), type I IFN signaling (e.g., ifnar1 and stat1), NF-ĸB-induced genes (e.g., tp53 and sod2), and ISGs (e.g., adar and isg20l2). It is important to note that the currently available X. laevis genome annotations remain incomplete, and immune genes are particularly under-annotated. Thus, while target predictions were performed with the most up-to-date annotation, our target prediction analysis may represent an underestimation of immune gene targets. Additionally, our analyses were restricted to 3′ UTRs, as although miRNAs have recently been shown to interact with 5′ UTRs and cds, the 3′ UTR is the location most often targeted by miRNAs (Lee et al. 2009; Fang and Rajewsky 2011). Nevertheless, targeting of these antiviral pathways is reflective of what has been observed in other vertebrates, where miRNAs that are differently expressed during immune responses have been found to target orthologs of these same immune genes (Boosani and Agrawal 2016). For example, miR-146a is induced by vesicular stomatitis virus in murine peritoneal macrophages (Hou et al. 2009) and targets irak2 in mice (Hou et al. 2009) and irf5 and stat1 in humans (Tang et al. 2009); miR-155 is induced during inflammatory responses in humans and mice (O’Connell et al. 2007; Tili et al. 2007) and targets ripk1 (Tili et al. 2007) and tab2 (Imaizumi et al. 2010); miR-208b and miR-499a-5p are induced in human hepatocytes in response to hepatitis C virus infection and target ifnar1 (Jarret et al. 2016); and miR-30a is induced in HeLa cells in response to coxsackievirus B3 and targets trim25 (Li et al. 2020). Identification of these differentially regulated miRNAs set the foundation for future studies to systematically validate the existence of physical functional interactions between differentially expressed X. laevis miRNAs and their putative antiviral target genes.

Our data sheds light on conserved post-transcriptional responses to viral dsRNA and the potential for involvement of miRNAs that target a core set of conserved antiviral genes. For example, it has been shown in mammals that poly(I:C) induces Sod2 expression to combat reactive oxygen species (ROS), which are associated with viral infection (Molteni et al. 2014). However, the mechanism through which poly(I:C) induces Sod2 remains unknown. We identified a miRNA (xla-miR-106-5p) that targets sod2 and is downregulated in response to poly(I:C). Similarly, human hsa-miR-106a-5p is also predicted to target Sod2. Thus, repression of xla-miR-106-5p may be a conserved mechanism through which Sod2 is induced to combat ROS accumulation in response to viral infection. In line with this, oxidative stress has been shown to repress miR-106-5p in humans (Tai et al. 2020), and downregulation of miR-106-5p has separately been shown to restore Sod activity to inhibit apoptosis and oxidative stress in rats (Li et al. 2017). Another interesting candidate is xla-miR-22-5p, which is upregulated by poly(I:C) treatment or FV3 infection and is predicted to target arf6 (ADP-ribosylation factor 6), a gene known to regulate TLR9-mediated immune signaling in response to exogenous DNA (Wu and Kuo 2012). As Arf6 enhances proinflammatory cytokine production (Wu and Kuo 2012), upregulation of xla-miR-22-5p and subsequent suppression of arf6 in response to viral infection may represent a strategy to modulate inflammation to prevent excessive cellular damage. Interestingly, arf6 is a conserved target of the miR-22-5p orthologs in humans, mice, rats, cows, and the Western clawed frog. Together, these observations suggest that miRNA-mediated mechanisms of regulating innate immunity are highly conserved, even in distantly related vertebrates.

FV3 is immunoevasive (Grayfer et al. 2015; Robert et al. 2017) and does not elicit robust antiviral responses in Xela DS2 skin epithelial-like cells (Bui-Marinos et al. 2021). Thus, we sought to profile changes in miRNA expression in response to FV3 for comparison to poly(I:C)-induced miRNA alterations to identify potential miRNA-based mechanisms of immunoevasion. Accordingly, FV3-infection is associated with changes in the expression of fewer miRNAs compared to poly(I:C), such that only 34 of the 65 miRNAs altered by poly(I:C) treatment were also affected by FV3 infection. Thus, the 31 miRNAs that are altered by poly(I:C) treatment and not FV3 infection, along with the 15 miRNAs that are altered by FV3 infection and not poly(I:C) treatment, may play a role in FV3’s ability to evade host immune systems in skin epithelial cells. In addition, some miRNA changes evident at 24 h post-poly(I:C) treatment were not evident until 72 h post-FV3 infection, and fewer of the differentially expressed miRNAs in FV3-infected Xela DS2 cells are predicted to target immune genes compared to those altered by poly(I:C). These findings suggest that FV3 infection may be associated with a dampened and (or) delayed immune response, which is in accordance with the theory that FV3 possesses immunoevasion capabilities.

Comparing the antiviral targets of poly(I:C)-altered miRNAs and FV3-altered miRNAs reveals several notable differences. For example, poly(I:C) treatment suppresses the expression of several miRNAs that are predicted to target multiple components of TLR7, TLR8, and TLR9 pathways, a layer of regulation that appears to be absent in FV3-infected cells. As FV3 is a dsDNA virus and dsDNA is a TLR9 ligand, the possibility exists that FV3 may actively interfere with miRNA-mediated regulation of TLR9 signaling as a facet of immunoevasion. It is also interesting to note that compared to FV3 infection, poly(I:C) treatment affects the expression of more miRNAs that target NF-kB-regulated genes and ISGs. This suggests that FV3 may actively interfere with the regulation of antiviral genes that are typically induced in response to the activation of antiviral signaling pathways. Despite these intriguing observations, it is important to recognize the caveat that extracellular delivery of poly(I:C) can initiate antiviral pathways that may be distinct from the pathways stimulated by FV3 infection (e.g., poly(I:C) can activate endosomal TLR3 pathways while viral dsRNA produced during FV3 replicative cycles may activate cytosolic sensing pathways such as RIG-I/MDA5). Additionally, numerous other viral components that are not mimicked by poly(I:C) (e.g., viral proteins, viral dsDNA) can be recognized by different PRRs. Thus, while extracellular poly(I:C) is not a perfect comparator, our previous studies have established that extracellular introduction of poly(I:C) activates antiviral pathways that lead to an effective antiviral response which protects against FV3 in Xela DS2 cells (Bui-Marinos et al. 2021) and serves to provide a reference for which miRNAs are regulated during an effective antiranaviral immune response. Future studies should examine miRNA profiles in response to poly(I:C) that is transfected into the cell, which may more closely approximate cytosolic PRR sensing of viral dsRNA replicative products.

To our knowledge, only one other study has simultaneously profiled changes in miRNA expression in response to poly(I:C) and an immunoevasive virus (porcine reproductive and respiratory syndrome virus (PRRSV)) (Wu et al. 2019). Similar to our observations, PRRSV infection of porcine alveolar macrophages resulted in alterations to only 33 of the 197 miRNAs found to be altered by poly(I:C) treatment (Wu et al. 2019). However, in contrast to the findings of the aforementioned study which did not detect any miRNAs altered by PRRSV infection that were not altered by poly(I:C) treatment, we detected differential expression of a handful of miRNAs in response to FV3 infection that were unaffected by poly(I:C) treatment. Of particular interest are two novel X. laevis miRNAs with no known orthologs in other species (novel-xla-miR-16-5p and novel-xla-miR-16-3p) that are strongly induced (150-fold and 342-fold, respectively, as measured by RT-qPCR) in response to FV3 but not poly(I:C). These novel miRNAs originate from the same pre-miRNA and are both predicted to target important genes involved in innate antiviral immunity, including tgfbr2 (transforming growth factor beta receptor 2; novel-xla-miR-16-5p) and trim14 (tripartite-containing motif 14; novel-xla-miR-16-3p). Tgfbr2 is a receptor for the well-studied TGF-β anti-inflammatory cytokine and viral induction of novel-xla-miR-16-5p may result in dampening of TGF-β signaling, preventing the host from limiting inflammatory responses to protect against cellular damage caused by excessive inflammation. Trim14 indirectly promotes type I IFN production by stabilizing the cytosolic dsDNA sensor cGAS (Chen et al. 2016). In the absence of Trim14, type I IFN antiviral responses are impaired (Chen et al. 2016). It is interesting to note that novel-xla-miR-16-3p was undetectable by RNA-seq in Xela DS2 cells that were not infected with FV3, but it is strongly induced in response to FV3 infection. As novel-xla-miR-16-3p targets a stabilizer of cGAS, which is a sensor of viral dsDNA, FV3 may force novel-xla-miR-16-3p expression to actively suppress cGAS signaling and evade recognition by the host immune system. Indeed, viral inhibition of cGAS signaling has been observed in other pathosystems (Eaglesham et al. 2019). Whether cGAS signaling is involved in responding to FV3 infection and the roles FV3-induced host miRNAs play in blocking cGAS signaling require further investigation.

Our study is the first to profile changes in host miRNA expression in response to ranavirus infection in skin cells modeled by a skin epithelial-like cell line, as the only other study that has profiled changes in host miRNAs in response to ranavirus infection was performed in spleen. Meng et al. (2018) profiled Andrias davidianus (Chinese giant salamander) miRNA responses to Chinese giant salamander iridovirus infection of spleen tissue and discovered numerous differentially expressed miRNAs predicted to target genes involved in important antiviral signaling pathways such as TLR, RIG-I/MDA5, and JAK-STAT pathways (Meng et al. 2018). Despite the similarities to pathways targeted in our study, only two miRNA families found to be differentially expressed in the study by Meng et al. (2018) were also found to be differentially expressed in response to FV3 in our study (miR-200 and miR-203, both downregulated). Thus, while the antiviral pathways targeted by miRNAs in response to different ranaviruses may mirror one another, it appears as though regulation of these pathways may be achieved using different miRNAs. Nevertheless, it will be important for future studies to verify that the differentially expressed miRNAs identified in our study function in antiviral responses and impact viral replication. As our study utilized a skin epithelial-like cell line, a potential limitation is that the findings may not be entirely representative of the complex processes occurring in whole skin tissue. Thus, to broaden our understanding of miRNA-mediated antiviral responses in frog skin, future studies are required to identify the contribution of frog miRNAs to antiviral responses in other cell types present in frog skin, as well as in whole skin tissue.

In addition to targeting host antiviral genes, several differentially expressed miRNAs are predicted to target viral transcripts corresponding to a myristoylated membrane protein (FV3 ORF 53 R), LCDV1 orf2-like protein, and hypothetical proteins. FV3 ORF 53R has been suggested to interact with the host cellular membrane and is required for virion assembly (Whitley et al. 2010), and thus may represent an ideal molecule for the host immune system to target. It is important to note that our target prediction analyses were performed with FV3 cds, as the current FV3 genome annotation lacks defined 3′ UTRs. While viral cds are frequently targeted by host miRNAs (Wang et al. 2012; Zheng et al. 2013; Ho et al. 2016), our FV3 target predictions may represent a fraction of the total FV3 mRNAs targeted by host miRNAs.

Recent studies have shown that transient inhibition of miRNA biogenesis is critical during the activation of IFN responses in HeLa cells, and leads to reduced levels of specific miRNAs, many of which target innate immune genes (Witteveldt et al. 2018). Accordingly, our study revealed that the majority of differentially expressed miRNAs that target immune genes are downregulated rather than upregulated, and both poly(I:C) treatment and FV3 infection suppressed dicer1 transcript levels. In line with this, studies have detected suppression of dicer1 but not drosha in response to viral infection of human, nonhuman primate, and murine cells (vaccinia virus infection of HeLa, Vero, and 3T3 cells) (Grinberg et al. 2012). Therefore, transient Dicer suppression may constitute a conserved antiviral response. Interestingly, we observed downregulation of drosha transcript levels in response to FV3 but not poly(I:C), which raises the possibility that Drosha suppression may be linked to immunoevasion tactics employed by FV3. Accordingly, studies have demonstrated that during infection of human 293T cells with Sindbis virus or vesicular stomatitis virus, Drosha is shuttled to the cytoplasm and repurposed to cleave viral RNA to inhibit viral replication (Shapiro et al. 2014). Thus, it stands to reason that immunoevasive viruses may be equipped to restrict Drosha function to facilitate viral replication. A survey of the literature revealed that direct manipulation of Drosha expression by a virus has, to the best of our knowledge, not yet been reported. Thus, it is unclear whether Drosha is a common target of viral immunoevasion or whether it may be a target that is unique to FV3.

Our sequencing efforts expanded the number of currently annotated X. laevis miRNAs from 247 to 380. Of the 133 novel X. laevis miRNAs we identified in Xela DS2 cells, 108 have never been detected in any species before. While several novel miRNAs were differentially expressed in response to poly(I:C) or FV3, the majority were not. A large portion of the novel miRNAs that were not altered by poly(I:C) or FV3 originate from repetitive regions found across multiple chromosomes and may thus function in controlling repetitive elements. Indeed, miRNAs have been recently characterized as emerging regulators of transposable elements (Pedersen and Zisoulis 2016). Thus, our findings are likely to stimulate exciting studies in research areas outside of innate immunity.

While our study significantly increased the number of annotated X. laevis miRNAs, it is important to note that our novel miRNA discovery efforts have likely not completed the map of the Xenopus miRNAome, as miRNAs are often exclusively expressed in specific cell types, tissue types, developmental stages, or cellular contexts (e.g., infection) (Bartel 2004). Additionally, our novel miRNA discovery pipeline involved removing reads with <10 counts in a given library, which may exclude the discovery of novel miRNAs expressed at low levels. However, removing reads with counts <10 improves the signal-to-noise ratio (Law et al. 2016) which enables the discovery of high confidence miRNAs (Friedländer et al. 2012), and this benefit arguably outweighs the limitation of missing some novel miRNAs.

To our knowledge, this study is the first to highlight the role of miRNAs in antiviral defenses within frog skin epithelial cells and sheds light on evolutionarily conserved post-transcriptional mechanisms of regulating antiviral responses in vertebrates. We identified poly(I:C)-induced antiviral responses involving miRNAs in frogs and uncovered differences between these responses and miRNA-mediated responses to the immunoevasive FV3. Our findings further our understanding of frog-FV3 interactions by identifying miRNAs as an added layer of antiviral gene regulation in frogs and uncovering the possibility that ranaviruses may hijack specific host miRNAs to facilitate immunoevasion. Future functional studies will be instrumental in solidifying these initial observations, and similar studies in other cell types and tissues will be crucial to uncovering cell type- and tissue-specific antiviral miRNA responses in frogs.