Temperature-dependent Small RNA Expression Depends on Wild Genetic Backgrounds of Caenorhabditis briggsae

Abstract Geographically distinct populations can adapt to the temperature conditions of their local environment, leading to temperature-dependent fitness differences between populations. Consistent with local adaptation, phylogeographically distinct Caenorhabditis briggsae nematodes show distinct fitness responses to temperature. The genetic mechanisms underlying local adaptation, however, remain unresolved. To investigate the potential role of small noncoding RNAs in genotype-specific responses to temperature, we quantified small RNA expression using high-throughput sequencing of C. briggsae nematodes from tropical and temperate strain genotypes reared under three temperature conditions (14 °C, 20 °C, and 30 C). Strains representing both tropical and temperate regions showed significantly lower expression of PIWI-interacting RNAs (piRNAs) at high temperatures, primarily mapping to a large ∼7 Mb long piRNA cluster on chromosome IV. We also documented decreased expression of 22G-RNAs antisense to protein-coding genes and other genomic features at high rearing temperatures for the thermally-intolerant temperate strain genotype, but not for the tropical strain genotype. Reduced 22G-RNA expression was widespread along chromosomes and among feature types, indicative of a genome-wide response. Targets of the EGO-1/CSR-1 22G-RNA pathway were most strongly impacted compared with other 22G-RNA pathways, implicating the CSR-1 Argonaute and its RNA-dependent RNA polymerase EGO-1 in the genotype-dependent modulation of C. briggsae 22G-RNAs under chronic thermal stress. Our work suggests that gene regulation via small RNAs may be an important contributor to the evolution of local adaptations.


Introduction
Temperature is a universal property of all ecosystems that influences diverse aspects of organismal fitness. In species that inhabit distinct temperature regimes across their geographic range, natural selection can lead different populations to evolve adaptations to the temperature conditions of their local habitat, creating within-species variation in the response to temperature stress. For example, populations of Drosophila melanogaster from tropical regions are more heat-tolerant than populations from cooler temperate regions (Rohmer et al. 2004). Similarly, populations of Arabidopsis thaliana found at high latitudes are more cold-tolerant than populations from low latitudes (Zhen and Ungerer 2008). Understanding the process of adaptation is a central theme of evolutionary biology, and so there is much interest in identifying the genetic basis of such local adaptations in diverse taxa (Ahrens et al. 2018). However, it is not clear in all cases what the genetic variants and molecular mechanisms underlying these local adaptations to temperature are, or how this variation contributes to differential fitness.
An excellent system for interrogating the molecular mechanisms of local adaptation is the nematode roundworm Caenorhabditis briggsae, a self-fertilizing hermaphrodite and close relative of the model organism Caenorhabditis elegans. Like C. elegans, C. briggsae has a widespread global range, and therefore has the potential to demonstrate local adaptation to varying habitats (Kiontke et al. 2011;Cutter 2015). In support of this possibility, C. briggsae isolates found at tropical and temperate latitudes form genetically distinguishable phylogeographic groups (Cutter et al. 2006;Félix et al. 2013). Examining how populations of C. briggsae found at different latitudes may show local adaptation, Prasad et al. (2011) demonstrated that strains of C. briggsae from three phylogeographic groups reared at varying temperatures showed distinct reproductive fitness responses. Strains from tropical latitudes had higher lifetime fecundity than temperate strains at high temperatures (30 °C) but lower fecundity than temperate strains at low temperatures (14 °C), consistent with local adaptation. Subsequent studies have further characterized how C. briggsae strains differ in their response to temperature, such as how strains differ genetically in behavioral responses to environmental temperature, including thermal preferences, isothermal tracking, motility, sperm and oocyte fertility, and acute cold sensitivity (Stegeman MBE et al. 2013;Poullet et al. 2015;Stegeman et al. 2019;Wang et al. 2021). mRNA sequencing of a tropical and a temperate strain reared at extreme temperatures found a significant genotype-temperature interaction for 56% of the differentially expressed genes (i.e., strains differed in temperature-dependent expression changes), many of which were associated with known biological processes such as translation and oogenesis (Mark et al. 2019). Many genes are upregulated in response to extreme cold stress in temperate, but not tropical, C. briggsae strains, and the role of a few of these genes in mediating strainspecific cold tolerance has been experimentally confirmed (Wang et al. 2021). Despite this work, the genetic mechanisms underlying these local adaptations in C. briggsae remain incompletely understood.
One potential contributor to locally-adaptive evolutionary responses that remains largely uncharacterized are small noncoding RNA regulatory pathways. Small RNAs (noncoding RNAs generally between 20 and 30 bases in length) associate with Argonaute proteins to regulate the expression of target loci, such as protein-coding genes and transposable elements, and include microRNAs, PIWI-interacting RNAs (piRNAs), and endogenous small interfering RNAs (endo-siRNAs) (Hoogstrate et al. 2014). Small RNA pathways can show adaptive differences between species (Palmer et al. 2018;Franchini et al. 2019), but the role of small RNA-mediated gene regulation in shaping differences between populations within a single species is less clear. However, there is some evidence linking small RNAs to local adaptations. Comparisons between different populations have revealed microRNA loci in nematodes (Jovelin and Cutter 2011), humans (Torruella-Loran et al. 2016), and plants (Xie et al. 2017) that show strong allele frequency differentiation that is consistent with local adaptation, though the functional consequences of this differentiation are not fully understood. The genes that carry out small RNA regulation may also show populationspecific signatures of positive selection, as seen for genes involved in the D. melanogaster piRNA pathway (Simkin et al. 2013). Sequence changes to small RNA binding sites can be involved in local adaptations as well, as has been reported for microRNA binding sites in populations of humans (Li et al. 2012) and D. melanogaster (Catalán et al. 2016). In addition to DNA sequence differences, natural populations can show variation in the expression of piRNAs (Ellison and Cao 2020), microRNAs (Arif et al. 2013), and siRNAs (Zhai et al. 2008). This variable expression between populations can cause intraspecific variation in morphology (Arif et al. 2013) or epigenetic regulation (Zhai et al. 2008), but the adaptive significance of this variation is usually uncharacterized. However, a few studies have identified differential microRNA expression between natural populations with known local adaptations, such as to heat stress (Graham and Barreto 2019) or hydrogen sulphide toxicity (Kelley et al. 2021), providing a link between microRNA-mediated gene regulation and locally-adaptive phenotypes. While much is unknown about how endo-siRNAs may be involved in local adaptation compared with microRNAs, the regulatory role of some endo-siRNAs in mediating responses to environmental stress (Wu et al. 2020) suggests that environmental differences between populations could lead to adaptive evolution in endo-siRNA pathways.
Studies in C. elegans indicate that small RNAs may play an important role in how these nematodes respond to temperature stress. For example, loss of 26G-RNAs (endo-siRNAs that are 26 bases long and begin with a 5′ guanine nucleotide) induces spermatogenesis defects at elevated temperatures (Conine et al. 2010). Additionally, production of piRNAs is repressed in response to increased temperatures (Belicard et al. 2018), and disruption of the piRNA pathway in C. elegans leads to fertility defects at high temperatures (Wang and Reinke 2008). It has been shown that some C. briggsae microRNAs have a high level of polymorphism between strains, and this variation may have functional consequences (Jovelin and Cutter 2011). Thus, differences in small RNA sequence and expression between tropical and temperate strains have the potential to contribute to the temperature-dependent mRNA transcriptome profiles and fitness differences observed in C. briggsae, although to date this idea has not yet been tested.
To explore this possibility, we characterize differential small RNA expression in response to chronic temperature stress for distinct temperate and tropical genotypes of C. briggsae. We observe genome-wide reductions in small RNA expression under high temperature stress that are specific to the temperate genotype, consistent with the reduced fitness of this strain at high temperatures. By classifying genes that are regulated by different small RNA pathways, we test how specific pathways may be responsible for this genotype-specific response. The temperaturedependent responses of distinct genotypes that we observe point to a possible role for the regulation of gene expression via small RNAs in the process of local adaptation.

Caenorhabditis briggsae Genotypes Show Contrasting Patterns of microRNA Expression at Low Temperatures
We investigated the effect of chronic temperature stress from rearing at 14 °C, 20 °C, and 30 °C on the small RNAs produced by different C. briggsae genotypes, using two strains that represent tropical (AF16) and temperate (HK104) latitude phylogeographic groups of this species. Endo-siRNAs with a 5′ guanine and 22-nucleotide length (22G-RNAs) make up a large fraction of the regulatory small RNAs in these samples ( fig. 1), which is consistent with other Caenorhabditis nematodes (Hoogstrate et al. 2014). On average, 31.4% of small RNAs corresponded to 22G-RNAs, accompanied by major a contribution of microRNAs (31.1%), as well as minor contributions of 21U-piRNAs (3.7%) and 26G-RNAs (1.4%) ( fig. 1,  supplementary material fig. S1, Supplementary Material MBE online). We found that 24.1% of all small RNAs aligned antisense to annotated protein-coding genes, with only 4.9% aligning in the sense orientation (supplementary material fig. S1, Supplementary Material online), as expected for endo-siRNAs such as 22G-RNAs. 22G-RNAs primarily aligned antisense to protein-coding genes (42.4% of all 22G-RNAs), pseudogenes (7.0%), repeats (5.5%), and transposable elements (4.1%), as well as unannotated regions of the genome (36.1%). About one quarter (26.1%) of all small RNAs aligned to regions of the genome without any functional annotation, though this is likely due to unannotated loci not present in our C. briggsae genome annotation. Small RNAs with a length of 24 nucleotides were strongly enriched for those starting with a 5′ cytosine (C) ( fig. 1), as has been seen previously in C. elegans and C. briggsae  Mean expression of small RNAs of different lengths that mapped perfectly to the Caenorhabditis briggsae genome for each genotypetemperature combination, separated by the 5′ nucleotide of mapped small RNAs. Pie charts give the distribution of genomic features that small RNAs mapped to. 'AS' indicates reads mapping antisense to a feature. MBE 2021). This enrichment was due to a single highly-expressed microRNA, cbr-mir-52 (WBGene00250149), located on chromosome IV. The strongest changes in microRNA expression were observed in temperate HK104 worms reared at 14 °C, which showed a 1.7-fold reduction in overall microRNA expression relative to HK104 worms at 20 °C (supplementary material fig. S1, Supplementary Material online). Interestingly, the tropical AF16 genotype showed the opposite pattern, with overall microRNA expression being 1.3-fold higher at 14 °C relative to 20 °C  (Shi et al. 2013;Beltran et al. 2019). We found that, overall, 21U-piRNAs show lower expression in worms reared at high temperature (30 °C) for both the tropical (AF16) and temperate (HK104) genotypes (figs. 1 and 2). This trend was driven almost entirely, however, by 21U-piRNAs that derived from the large piRNA gene cluster on the left arm of chromosome IV High Temperatures Reduce 22G-RNA Expression Only in the Temperate Genotype In addition to identifying some shared temperaturedependent expression responses, we also observed genotype-specific small RNA expression changes in response to temperature. Using multidimensional scaling, expression profiles of 22G-RNAs that mapped antisense to annotated genomic features clustered into three distinct sample types: AF16 animals from all temperatures, HK104 worms from 14 °C and 20 °C, and a separate group of HK104 reared at 30 °C (supplementary material fig. S3, Supplementary Material online). Using linear modeling to compare expression between genotype-by-temperature combinations, temperate HK104 animals raised at 30 °C showed a 2.3-fold decrease in 22G-RNA reads aligning antisense to protein-coding genes relative to HK104 animals raised at 20 °C (P = 0.014), despite tropical AF16 animals having on average nearly identical gene-antisense 22G-RNA expression between 20 °C and 30 °C (P = 0.956; fig. 1 and fig. 3, supplementary material fig. S1, Supplementary Material online). This represented a significant genotype-temperature interaction on 22G-RNA expression at 30 °C (P = 0.050), consistent with especially strong perturbation of 22G-RNA expression in HK104 at high temperatures. This trend also held for 22G-RNAs aligning antisense to pseudogenes (2.3-fold lower in HK104 at 30 °C; P = 0.016) ( fig. 3), and was nominally similar if non-significant for repeat regions (1.6-fold lower; P = 0.338) and transposable elements (1.4-fold lower; P = 0.345), indicative of a genome-wide rather than locusspecific reduction in 22G-RNA expression for the temperate HK104 genotype reared at 30 °C.
Supporting the idea of a global regulatory difference, we observed that expression of all 22G-RNAs across all six C. briggsae chromosomes decreased 1.5-fold for HK104 animals reared at 30 °C compared with 20 °C, with decreases ranging from 1.3-fold on chromosome III to 2.0-fold on chromosome II ( fig. 4, supplementary material fig. S4, Supplementary Material online). Linear modeling of the 22G-RNAs aligning to each chromosome again found that, accounting for differences between chromosomes, this genome-wide decrease at 30 °C is statistically significant (P = 0.0013). Tropical AF16 animals, by contrast, showed 1.2-fold higher average 22G-RNA expression at 30 °C than at 20 °C on all chromosomes (ranging from 1.03-fold higher on chromosome I to 1.3-fold higher on chromosome IV) except chromosome II, which had 1.1-fold lower average expression ( fig. 4), although this overall increase at 30 °C was not statistically significant (P = 0.161). These contrasting patterns between genotypes again resulted in a significant genotype-temperature interaction at 30 °C (P = 9.3 × 10 −4 ). 22G-RNA expression was significantly higher on the X chromosome compared with all other chromosomes, regardless of genotype or temperature treatment (P < 6.1 × 10 −8 ; fig. 4).
Dividing the genome into bins 100 kb wide, 82-90% of bins across HK104 chromosomes showed this pattern of depressed 22G-RNA expression at 30 °C versus just 52-78% of bins for AF16 (table 1). Moreover, the mean magnitude of decreased expression was larger for HK104 than AF16 on all six chromosomes (table 1). Both the chromosome arm and center domains reflected these patterns (supplementary material tables S1-S2, Supplementary Material online), although centers showed a more pronounced difference between genotypes: at 30 ° C, 35.4% of bins on chromosome centers had lower 22G-RNA expression in HK104 but not AF16, compared with only 19.5% of bins on chromosome arms (Fisher's exact test, P = 9.6 × 10 −9 ).
Using the 22G-RNA expression profiles of 14,283 genomic features to which 22G-RNA sequence reads aligned with antisense orientation (including protein-coding MBE genes, pseudogenes, repeats, and transposable elements), we detected significant differential expression for 8,531 of them due to strain genotype, rearing temperature, or a non-additive genotype-temperature interaction (i.e., expression changes due to temperature that differed between the AF16 and HK104 genotypes) (supplementary material fig. S5  Given the post-transcriptional regulatory role of 22G-RNAs on mRNA for the genes from which they both originate, we examined the association between 22G-RNA and mRNA expression levels across genotype

MBE
and temperature treatment combinations. We could directly relate 22G-RNA and mRNA expression for 12,623 genes with data that derived from the same total RNA samples [22G-RNA from this study and mRNA from Mark et al. (2019)]. Overall, these 12,623 genes show moderate but significant positive correlations between 22G-RNA and mRNA expression across all genotype-temperature combinations (P < 2. Upon testing for a shared profile of genotype only, temperature only, or genotype-by-temperature interaction for 22G-RNA and mRNA expression of a given protein-coding gene, we detected surprisingly little overlap (Supplementary material fig. S9, Supplementary Material online). For example, only 30.9% (1269/4108) of genes (with both 22G-RNA and mRNA expression data available) showing a 22G-RNA genotype-temperature interaction also showed an mRNA genotype-temperature interaction. This overlap was marginally statistically insignificant (hypergeometric test, P = 0.055), and the overlap for other differential expression categories was even lower (supplementary material fig. S9, Supplementary Material online). These observations point to the possibility that additional regulatory controls over 22G-RNA and mRNA FIG. 3. Total expression of 22G-RNAs aligning antisense to protein-coding genes, pseudogenes, repeat regions, and transposons for each genotype-temperature combination, considering only the genes, pseudogenes, repeat classes, and transposon classes that are present in the genome annotations for both strains. Expression is in units of RPM. Error bars give ±1 standard deviation around the mean of each set of replicates.

MBE
expression complicates their genotype-by-temperature interactions as reflected in RNA abundances.
Reduced 22G-RNA Expression is Specific to the EGO-1/CSR-1 Pathway RNA-dependent RNA polymerases (RdRPs), such as EGO-1 and RRF-1, initiate the biogenesis of 22G-RNAs that then bind to Argonaute proteins, including CSR-1, WAGO-1, and HRDE-1. These Argonautes serve as effectors that modulate expression of target genes in the germline, mediated by base complementary binding with 22G-RNA sequences. CSR-1 binds 22G-RNAs made by EGO-1, WAGO-1 binds 22G-RNAs made predominantly by RRF-1 and by EGO-1 to a lesser extent, and HRDE-1 binds To investigate these 22G-RNA pathway proteins as potential sources of temperature-dependent reductions in 22G-RNA expression, we quantified their mRNA expression profiles. Cbr-ego-1, Cbr-rrf-1, Cbr-csr-1, Cbr-wago-1, and Cbr-hrde-1 all shared similar mRNA expression across all six genotype-temperature combinations ( fig. 5A). In the tropical AF16 strain, mRNA expression of these five genes was insensitive to rearing temperature and, at 20 °C, mRNA expression was indistinguishable between AF16 and HK104. The stressful rearing temperatures of both 14 °C and 30 °C, by contrast, typically led to reduced expression of these RdRPs and Argonautes in the temperate HK104 strain ( fig. 5A). We also profiled the mRNA expression of several genes in the 26G-RNA and piRNA pathways, which act upstream of 22G-RNA production (Hoogstrate et al. 2014), although these did not generally show a  These observations point to the possibility that the differences in 22G-RNA expression (figs. 3-4, Supplementary Material online) arise from changes in mRNA expression of genes in the 22G-RNA pathway. Reductions to 22G-RNA expression should be strongest in the genes targeted by a given pathway if changes to the pathway's Argonaute or RdRP expression cause lower 22G-RNA expression. To explore this possibility, we first compiled sets of genes targeted by the 22G-RNAs associated with EGO-1, RRF-1, CSR-1, WAGO-1, or HRDE-1 from previous RNA-sequencing studies in C. briggsae (for CSR-1) or C. briggsae orthologs of C. elegans targets (for EGO-1, RRF-1, WAGO-1, HRDE-1, and the CSR-1a and CSR-1b isoforms of CSR-1) Vasale et al. 2010;Buckley et al. 2012;Tu et al. 2015;Charlesworth et al. 2021). The resulting putative target sets contain between 146 and 4,839 genes (median 1,706 genes), for which we obtained C. briggsae 22G-RNA expression data of 121-3953 genes (median 1,456 genes) (table 2; supplementary  material table S4, Supplementary Material online). Consistent with their function in shared 22G-RNA pathways, EGO-1 and CSR-1 show extensive target overlap (1,427 out of 1,706 EGO-1 targets are also CSR-1 targets, 83.6%), as do RRF-1 and WAGO-1 (72 out of 146 RRF-1  fig. S12, Supplementary Material online), consistent with the known roles of CSR-1 and WAGO-1 in promoting and silencing gene expression, respectively, via 22G-RNA targeting Conine et al. 2013).
Comparing the 22G-RNA expression of the targets of EGO-1/CSR-1, HRDE-1, and RRF-1/WAGO-1 between genotypes and temperatures revealed striking differences between the pathways. Genes targeted by EGO-1 and CSR-1 showed 22G-RNA profiles that mirrored the genome-wide pattern: compared with 20 °C, 22G-RNA expression was significantly reduced specifically at 30 °C in HK104 but not AF16 worms ( fig. 6, Supplementary Material online), with AF16 worms showing slightly increased expression at 30 °C. While both CSR-1a and CSR-1b targets shared this same trend, it was most pronounced for CSR-1b targets, and genes exclusively targeted by CSR-1a and not CSR-1b did not show this trend at all (supplementary material fig. S13, Supplementary Material online). These observations implicate the germline-constitutive CSR-1b short isoform, and not the intestinal and spermatogenesis-associated CSR-1a long isoform (Charlesworth et al. 2021), in this pattern of genotype-dependent differential expression of 22G-RNAs at 30 °C.
Genes targeted by HRDE-1 showed the same overall pattern as EGO-1 and CSR-1 targets, but weaker in magnitude, with high temperatures leading to significantly reduced 22G-RNA expression in HK104, but not AF16 ( fig. 6). This pattern is likely due to partial overlap between EGO-1 and HRDE-1 targeting, as 37.5% (354/945) of HRDE-1 targets are also EGO-1 targets. Indeed, partitioning HRDE-1 targets into those genes also targeted by EGO-1 and those not targeted by EGO-1 showed that the decrease in HK104 worms at 30 °C was strongest in HRDE-1 targets that are also EGO-1 targets (supplementary material fig. S14, Supplementary Material online). Surprisingly, genes targeted by RRF-1 and WAGO-1 did not show the genome-wide pattern. Instead, 22G-RNA expression was either the same (RRF-1 targets) or higher (WAGO-1 targets) at 30 °C in HK104 compared with AF16 ( fig. 6), though at 30 °C RRF-1 targets showed a significant decrease in HK104 but not AF16, similar to HRDE-1 targets. These profiles of target gene 22G-RNA expression match the profiles of target gene mRNA expression across genotypes and temperatures, including a broad decrease in mRNA expression for EGO-1 and CSR-1 targets at 30 °C in HK104 (supplementary material fig. S15, Supplementary Material online), consistent with the overall positive association between 22G-RNA and mRNA expression levels. Because temperature affects developmental timing, and gene expression varies across development, it is conceivable that differences in 22G-RNA expression between genotypes and temperatures might arise from subtle differences in developmental stage among young adults. To address this possibility, we examined the 22G-RNA expression for orthologs of RdRP and Argonaute target genes that show non-dynamic gene expression across development in C. elegans (Hendriks et al. 2014). This subset of genes still showed significantly decreased 22G-RNA expression of EGO-1/CSR-1 targets at 30 °C in HK104 but not AF16 (supplementary material fig. S16, Supplementary Material online), suggesting that this pattern is unlikely to be due to differences in developmental staging.
Our gene-level counts of 22G-RNA expression were normalized without using the standard trimmed mean of M-values (TMM) method to avoid the assumption that most genes do not show differential expression (Evans et al. 2018), as our data violated this assumption. Nevertheless, these results are robust to TMM normalization. Applying TMM scaling to the 22G-RNA expression of these 14,283 genes, we find that targets of EGO-1 and CSR-1 (but not HRDE-1) still show greatly decreased 22G-RNA expression in HK104 but not AF16 worms at 30 °C despite the normalization removing this pattern from the background distribution of genes (supplementary material fig. S17, Supplementary Material online). Together, these results support the idea that temperature stress at 30 °C uniquely disrupts the germline-expressed EGO-1 RdRP/CSR-1 Argonaute pathway in the temperate HK104 genotype. This occurs primarily through the CSR-1b isoform, and depresses the expression of hermaphrodite germline-associated 22G-RNAs and their target genes.

Discussion
Strains of C. briggsae from tropical and temperate latitudes show responses to temperature stress that are consistent with local adaptation: tropical strains are more fecund than temperate strains when raised at high temperatures, MBE whereas temperate strains are more fecund and viable than tropical strains when exposed to low temperatures (Prasad et al. 2011;Wang et al. 2021). Their responses are also genetically distinct in terms of behavior (Stegeman et al. 2013;Stegeman et al. 2019), mRNA expression (Mark et al. 2019;Wang et al. 2021), and in sperm fertility and germline development (Prasad et al. 2011;Poullet et al. 2015). Here, we incorporate differentiation in the expression of endogenous small RNAs by sequencing small RNAs from tropical and temperate C. briggsae genotypes raised under distinct temperature conditions. We find that some temperature-induced changes in small RNA expression are common to both genotypes, such as reduced piRNA expression, whereas other changes are genotype-specific. In particular, we document how globally reduced 22G-RNA expression is specific to the temperate HK104 genotype at a high 30 °C rearing temperature, possibly helping to explain why high temperatures have deleterious fitness consequences on this genetic background. Thus, we find that the evolution of small RNA pathways may be a contributor to the local adaptations to temperature seen in C. briggsae, an idea that motivates experiments to test explicitly for causality. While much is unknown about the role of small RNAs (and endo-siRNAs in particular) in local adaptation, our work suggests that these pathways are potentially important and warrant further study in this and other species.
Both the tropical AF16 genotype and the temperate HK104 genotype showed decreased overall piRNA expression at high temperatures, in line with similar findings in C. elegans and other C. briggsae strains (Belicard et al. 2018). The three distinct genomic clusters of piRNA loci (Shi et al. 2013;Beltran et al. 2019), however, showed differing expression responses to extreme temperatures. The piRNAs encoded in the cluster on the left arm of chromosome IV (0.6-7.0 Mbp) showed decreased expression at 30 °C, whereas those on the right arm (13.2-15.0 Mbp) showed increased expression at 14 °C. It is currently not known whether these clusters participate in different silencing pathways or target different sets of transcripts in C. briggsae, and so further work is needed to determine the full significance of this cluster-dependent piRNA expression. Cluster-specific differences in piRNA expression at varying temperatures could theoretically be caused by differences in chromatin organization, as piRNA loci in C. elegans normally are highly associated with regions of repressive chromatin (Beltran et al. 2019), though the chromatin organization of these clusters in C. briggsae at different temperatures is currently unknown. In nematodes, piRNAs act in the germline to silence transposons and other foreign DNA sequences (Lee et al. 2012), and disruptions of this pathway cause fertility defects at high temperatures (Wang and Reinke 2008). However, the role of piRNAs in cold temperature FIG. 6. 22G-RNA expression of genes targeted by EGO-1, RRF-1, CSR-1, WAGO-1, or HRDE-1, as well as all genes showing 22G-RNA expression, for each genotype-temperature combination. Significant differences between temperature treatments of the same genotype were determined using a Wilcoxon rank sum test, with *P < 0.05, **P < 0.01, ***P < 0.001, and n.s. indicating non-significance.

MBE
phenotypic responses has yet to be characterized. The cluster-specific expression changes of C. briggsae piRNAs at extreme temperatures in both genotypes may therefore point to a role associated with the general reduction in lifetime fecundity seen at these temperatures (Prasad et al. 2011).
For 22G endo-siRNAs, the most common type of small RNA in our samples, we observed the largest expression changes for the temperate HK104 strain when reared at 30 °C. We documented substantial genome-wide reductions in 22G-RNA expression in this genotype-by-temperature combination for multiple different types of loci, on multiple different chromosomes, and at numerous different loci across the full length of each chromosome. Together, these observations point to a global, genome-wide decrease in expression that is specific to the HK104 genotype when reared under chronic high thermal stress. These broad changes in 22G-RNA expression may contribute to the significant decrease in lifetime fecundity seen in the temperate HK104 strain when reared at 30 °C (Prasad et al. 2011); disruptions to C. briggsae 22G-RNA pathways are known to have a negative effect on offspring viability (Tu et al. 2015), though the possible contribution of 22G-RNA misexpression to HK104's heat intolerance requires further functional study. Evidence from mRNA sequencing suggests that many genes are indeed misregulated in the HK104 strain under heat stress (Mark et al. 2019), although it remains unclear to what extent perturbed mRNA expression levels depend on transcriptional versus post-transcriptional regulatory mechanisms, such as 22G-RNA-associated Argonaute pathways.
On average, genes targeted by the Argonaute CSR-1 and its RNA-dependent RNA polymerase (RdRP) EGO-1 showed reduced expression of antisense 22G-RNAs at high temperatures in the temperate HK104 strain but not the tropical AF16 strain. We generally did not observe this genotype-temperature interaction, however, for targets of the RRF-1/WAGO-1 pathway. This difference in 22G-RNA target profiles may arise from differences in gene expression of the RdRPs: Cbr-ego-1 showed substantially reduced mRNA expression at high temperatures in only HK104, whereas Cbr-rrf-1 showed a weaker reduction. As EGO-1 is needed for the biogenesis of the 22G-RNAs bound by CSR-1 ), downregulation of EGO-1 expression in this genotype under heat stress should lead to the production of fewer 22G-RNAs that can bind to CSR-1. However, given that the temperate genotype at low temperatures also showed similarly reduced EGO-1 gene expression without also showing reduced EGO-1/CSR-1 22G-RNA expression, EGO-1 gene expression alone cannot explain the reduction of 22G-RNAs at high temperatures. It is possible that this reduction is instead driven by reduced CSR-1 gene expression (e.g., if 22G-RNAs that are not bound to CSR-1 are more likely to be degraded), as CSR-1 in the temperate strain showed a stronger decrease in gene expression at high temperatures than at low temperatures, though other components of the EGO-1/CSR-1 pathway not studied here may also be involved. We also observe decreased mRNA expression for EGO-1/CSR-1 targets in HK104 at high temperatures, as expected given that CSR-1 plays a positive role in promoting gene expression of its targets (Conine et al. 2013;Wedeles et al. 2013).
The overall positive association between 22G-RNA expression and mRNA expression, however, makes it difficult to conclude whether decreased 22G-RNA expression causes decreased mRNA expression or vice versa. For instance, at 30 °C, targets of WAGO-1 show increases in both 22G-RNA expression and mRNA expression in HK104 compared with AF16, despite WAGO-1's role in gene silencing ). Nevertheless, large-scale misregulation of EGO-1/CSR-1 targets in the germline may partially explain the reduced sperm fertility that temperate C. briggsae strains show at elevated temperatures (Prasad et al. 2011). If so, then temperate C. briggsae strains may have acquired mutations involving the EGO-1/CSR-1 pathway that compromise the adaptive high temperature tolerance of tropical strains and ancestral genotypes-perhaps by causing the downregulation of EGO-1 or CSR-1 at high temperatures-although further study is needed to test possible mechanisms directly.
C. briggsae and other Caenorhabditis nematodes express both a long isoform (CSR-1a) and a short isoform (CSR-1b) of CSR-1, which differ due to an additional 163 N-terminal amino acids in CSR-1a (Tu et al. 2015;Charlesworth et al. 2021). We found that HK104 worms showed reduced 22G-RNA expression at 30 °C compared with 20 °C in targets of both CSR-1 isoforms, though this trend was more pronounced in CSR-1b targets than CSR-1a targets. In C. elegans, genes targeted by only the CSR-1b isoform are enriched for EGO-1 targets whereas genes targeted by only the CSR-1a isoform are enriched for RRF-1 targets (Charlesworth et al. 2021;Nguyen and Philips 2021). Consistent with our conclusion that EGO-1 and not RRF-1 is the key RdRP associated with the genotypespecific 22G-RNA reduction at 30 °C, we did not observe this 22G-RNA reduction in the genes targeted only by CSR-1a. Interestingly, expression of CSR-1a is required for full fertility under high temperature stress in C. elegans due to its role in spermatogenesis (Charlesworth et al. 2021). As CSR-1a is only expressed in the germline during the L4 larval stage when spermatogenesis occurs (Charlesworth et al. 2021), it is possible that our RNA sequencing of adult worms has missed temperaturedependent differences between strains in the CSR-1a pathway associated with spermatogenesis.
One limitation to our analyses of the different 22G-RNA pathways is that the sets of EGO-1 and RRF-1 targets used here are almost certainly incomplete. For instance, 70.0% (3388/4839) of the genes we defined as CSR-1 targets were not also defined as either EGO-1 or RRF-1 targets, despite the fact that EGO-1 is thought to be the main source of 22G-RNAs bound by CSR-1 ). Part of this incompleteness is simply a byproduct of the experimental approaches used to define them, as the single-gene mutants used may not be sufficient in recovering the full set of EGO-1 or RRF-1 target genes, given MBE that EGO-1 and RRF-1 are partially redundant in 22G-RNA biogenesis Gu et al. 2009;Vasale et al. 2010). In addition, we relied on orthology of gene targets between C. elegans and C. briggsae for most Argonaute pathway members. Comprehensive profiling of EGO-1 and RRF-1 targets in C. briggsae itself will benefit future analysis for how these RdRPs may influence C. briggsae local adaptation.
We found that 22G-RNA expression correlated positively with mRNA expression in all genotype-temperature combinations, even among genes targeted by WAGO-1, which is thought to silence its targets ). This general positive correlation is likely due to the requirement of mRNAs for the biogenesis of complementary 22G-RNAs, as these mRNAs serve as a template utilized by RdRPs in 22G-RNA production (Hoogstrate et al. 2014). However, we found that genes with the highest levels of 22G-RNA expression, including WAGO-1 targets, negatively correlated with mRNA expression, consistent with 22G-RNA-mediated silencing at sufficient levels of 22G-RNA expression. This nonlinear relationship between 22G-RNA and mRNA expression also has been observed in the C. elegans germline (Bezler et al. 2019).
Surprisingly, however, we observed little overlap between genes with genotype-or temperature-dependence in 22G-RNA expression and those genes showing the same effects for mRNA expression. This disparity could suggest that most changes in 22G-RNA expression are not strong enough to perturb mRNA expression in a way that matches the 22G-RNA profile, especially if pretranscriptional regulation (e.g., chromatin modifications that regulate transcription initiation) is a potent factor in controlling genotype-or temperature-dependent changes in transcription. Another potential explanation is that 22G-RNAs may disproportionately regulate target gene expression through processes downstream of mRNA. For example, endo-siRNAs in plants that are 22 nucleotides in length are capable of repressing the mRNA translation of target genes (Wu et al. 2020). Similar translation inhibition may also occur in some genes showing 22G-RNA effects in C. briggsae, although data supporting a link between 22G-RNAs and translational regulation in Caenorhabditis are lacking.
Overall, we find that tropical and temperate strains of C. briggsae show significant differences in how endogenous small RNA pathways respond to the same chronic temperature stress. Heat stress induced a genome-wide decrease in the expression of endogenous 22G-RNAs, but only in the strain derived from a temperate latitude. We hypothesize that genotype-specific differences in 22G-RNA expression, potentially caused by altered expression of genes involved in CSR-1b-associated 22G-RNA biogenesis, may contribute to heightened sensitivity to heat stress with profound consequences for animal fitness. C. briggsae thus provides a powerful experimental system to elucidate how small RNA mechanisms may contribute to local adaptations in nature.

Small RNA Sequencing Datasets
To assess the sensitivity of small RNA expression to rearing temperature and genetic background, we isolated and sequenced small RNAs for two isogenic strains of C. briggsae, each reared under defined temperature conditions throughout development (14 °C, 20 °C, and 30 °C). The AF16 strain represents a tropical latitude genotype whereas HK104 represents a genetically differentiated temperate latitude genotype, with the extreme 14 °C and 30 °C temperatures representing chronic stressful conditions for both genotypes (Cutter et al. 2006;Prasad et al. 2011;Félix et al. 2013;Thomas et al. 2015). Worms were cultured as described in Mark et al. (2019). Briefly, worms were plated on 100 mm Nematode Growth Medium plates with 1.5 mL of concentrated OP50 Escherichia coli. Plates were incubated at 14 °C, 20 °C, and 30 °C for approximately 150, 65 and 48 h, respectively, to the mature adult stage. Total RNA from gravid adult hermaphrodites in mass isogenic cultures of each strain-temperature combination was isolated with TRIzol extraction and isopropanol precipitation (Tu et al. 2015) from three biological replicates, for a total of 18 distinct samples that correspond to the identical samples analyzed for mRNA expression by Mark et al. (2019). Using the Ambion MirVana kit, we separated the small RNA fractions <200 bp and then isolated small RNAs between 18 and 30 nucleotides using gel extraction from 15% acrylamide gels, as adapted from Claycomb et al. (2009) and Gu et al. (2011). The samples were treated with tobacco acid pyrophosphatase to convert triphosphorylated RNA into a monophosphorylated form prior to ligating the 3′ linker, followed by 15% polyacrylamide/7 M urea gel extraction Gu et al. 2011). A 5′ linker sequence, possessing one of nine 4-nucleotide barcodes at its 3′ end, was ligated to the samples to enable multiplexing. This was followed by reverse-transcription cDNA synthesis, and PCR amplification, to add Illumina sequencing adapters with an additional 8-nucleotide indexing barcode on the 3′ adapter.
Small RNA libraries from each replicate were sequenced at the Tufts University Core Facility. Samples containing 202,528,982 total raw reads from AF16 and HK104 were demultiplexed and had their 5′ 4-nucleotide barcodes removed using a custom Perl script. Before barcode removal, Cutadapt v1.2.1 (Martin 2011) was used to remove 3′ adapters (parameters: -e 0.1 -m 21 -M 34 -a CTGTAGGCACCATCAAT) and discard reads that contained at least 14 nucleotides of the 5′ adapter (parameters: -e 0.1 -untrimmed-output=<output file > -O 14 -g TCTACAGTCCGACGATC). Reads were discarded if, after adapter and barcode trimming, the read was longer than 30 nucleotides or shorter than 17 nucleotides. After trimming, the 18 samples averaged 9,660,349 total reads each, giving an average of 28,981,046 reads for each of the six treatments (supplementary material table S5, Supplementary Material online). Analysis of all samples using FastQC v0.11.8 (Andrews 2010) confirmed for 17 MBE of them that FASTQ files contained, as expected, a majority of reads of approximately 22 bases in length and no detected adapter sequences. However, one replicate of HK104 grown at 30 °C showed an unusually low read count (1,298,689 vs. >4,500,000 reads in all other samples; supplementary material table S5, Supplementary Material online) with a length distribution skewed heavily towards short 17-19 nucleotide reads. These anomalies led us to exclude this sample from downstream analyses other than genome alignment (see below), resulting in an average of 10,152,211 reads per sample and 28,764,597 reads per treatment.
Generation of a Genotype-specific Reference for HK104 Worms To align these small RNAs to the C. briggsae genome for expression quantification, we used the reference genome sequence and gene annotations from WormBase release WS272 (Lee et al. 2018). Locations of annotated C. briggsae piRNA loci were taken from the file 'cbg_piRNAs_motifs_and_21Us.bed' from supplementary data S1, Supplementary Material online of Beltran et al. (2019). Since this reference genome is based on the tropical AF16 strain, we also generated a genotype-specific reference for the temperate HK104 strain. Whole-genome sequencing reads for the HK104 strain from Mark et al. (2019) were obtained from NCBI's Sequence Read Archive (accession SRR8333803) and aligned to the AF16 reference genome. Alignments were performed with NextGenMap v0.5.2 (Sedlazeck et al. 2013) in paired-end mode, with all other parameters left as default, including the default cutoff of 65% alignment identity to consider a read as successfully mapped. NextGenMap is able to account for a high level of divergence between the sequences being aligned, making it suitable for performing alignments between different strains. These alignments were used to call single nucleotide polymorphisms (SNPs) and short indels between the AF16 and HK104 genotypes with the HaplotypeCaller function from the GATK v4.1.4.1 software toolkit run with default parameters (McKenna et al. 2010). Using a custom Python script, we substituted these SNPs and indels into the AF16 reference genome sequence in order to create an HK104 'pseudo-reference' genome sequence. We then ran flo (Pracana et al. 2017) with default parameters to create a gene annotation file for HK104 by transferring gene annotations (microRNA genes, protein-coding genes, pseudogenes, transposable elements, repeat regions, rRNA genes, tRNA genes, piRNA genes, and other noncoding RNA genes)from the WormBase AF16 reference to the coordinates of the HK104 pseudo-reference sequence (available at https:// github.com/Cutterlab/C_briggsae_small_RNAs). HK104 protein-coding gene annotations were kept only for 17,038 genes (81.8% of 20,829 protein-coding genes total) in which all constituent parts (e.g., exons, UTRs) could be successfully converted from AF16 coordinates to HK104 coordinates. piRNA gene annotations were only kept if the resulting piRNA gene was still 21 bases long and started with a T nucleotide in the HK104 pseudo-reference.

Alignment of Reads to Reference Genomes
Small RNAs were aligned to the C. briggsae genome using the STAR v2.7.2a read alignment software (Dobin et al. 2013). All small RNA samples were aligned to both the AF16 genome (using the genome sequence and gene annotations from WormBase) and the HK104 genome (using the genotype-specific genome sequence and gene annotations created here). Both genomes were compiled with STAR's genomeGenerate function using the following options: -genomeSAindexNbases 12 -sjdbOverhang 29. Small RNA alignments to each genome were run with the following options: -quantMode GeneCounts -limitBAMsortRAM 30000000000 -outSAMtype BAM SortedByCoordinate -outFilterMultimapNmax 50 -outFilterMultimapScoreRan ge 0 -outFilterMismatchNoverLmax 0.05 -outFilterMa tchNmin 16 -outFilterScoreMinOverLread 0 -outFilte rMatchNminOverLread 0 -alignIntronMax 1. AF16 small RNA reads mapped best to the AF16 reference genome and HK104 reads mapped best to our HK104 pseudo-refe rence genome (supplementary material

Counting Aligned Small RNA Reads
Aligned small RNA reads from all samples were counted using a custom R pipeline (Charlesworth et al. 2021; https://github.com/ClaycombLab/Charlesworth_2020). Mapped reads were counted only if they aligned to the genome (either sense or antisense) without any insertions, deletions, or mismatches, but allowing for gaps due to introns. Each perfectly-aligned read was annotated with its length, first nucleotide, and coordinate location in the genome, using the location of the primary alignment (as determined by STAR) for reads that aligned to multiple locations. For each annotated genomic feature on C. briggsae's six chromosomes (microRNA genes, protein-coding genes, pseudogenes, transposable elements, repeat regions, rRNA genes, tRNA genes, piRNA genes, other noncoding RNA genes, and a single feature comprising all unannotated regions), we counted how many perfectly-aligned reads of each length and first nucleotide aligned completely to this feature in each replicate, separately for sense alignments and antisense alignments (supplementary material table S5, Supplementary Material online). For reads that aligned perfectly to multiple different locations in the genome at the same time, each alignment of that read to a feature was only counted as 1/N reads for that feature, where N is the total number of times the read aligned to a feature of that class (e.g., protein-coding genes).

MBE
These counts were used to determine the total number of microRNAs, 21U-piRNAs, and 22G and 26G endo-siRNAs in each replicate expressed from each C. briggsae chromosome, including unplaced genomic scaffolds. We defined microRNAs as reads aligning to the sense strand of annotated microRNA genes, 21U-piRNAs as reads with a length of exactly 21 bases beginning with a uridine (U) nucleotide aligning sense to annotated piRNA genes, and 22G-and 26G-RNAs as any reads with a length between 21-23 bases and 25-27 bases, respectively, and beginning with a guanine (G) nucleotide. Due to their observed numerical abundance, we primarily focused on characterizing the expression of 22G-RNAs. The total number of 22G-RNAs aligning to each of the six C. briggsae chromosomes (excluding unplaced genomic scaffolds) was fit using the linear model: 22G-RNA expression = chromosome of origin effect + genotype effect + temperature effect + genotype:temperature interaction effect. This linear model was used to assess the significance of the temperature effect for each genotype and the genotype-temperature interaction, at 30 °C relative to 20 °C, as well as the chromosome effect for each chromosome. All statistical analyses were done in R v4.1.1 (R Core Team 2021). Coordinates for chromosome arm and center domains in the current C. briggsae genome assembly (cb4) were taken from Ross et al. (2011), combining the short chromosome tip regions with the adjacent arm domains.

22G-RNA Expression Analysis
To investigate the expression of 22G endo-siRNAs across genotypes and temperature treatments, we quantified the number of 22G-RNAs aligning antisense to all annotated repeats, transposons, pseudogenes, and proteincoding genes in the C. briggsae genome, using only the 16,916 unique features whose coordinates completely transferred from the AF16 genome annotation to the HK104 genome annotation. Repeat and transposon classes were kept if at least one copy of the repeat/transposon was present in both genome annotations (i.e., without requiring every copy to be present in both genomes). The total number of 22G-RNAs aligning antisense to these 16,916 features was fit, separately for protein-coding genes, pseudogenes, repeats, and transposons, using the linear model: 22G-RNA expression = genotype effect + temperature effect + genotype:temperature interaction effect. These linear models were used to assess the significance of the temperature effect for each genotype and the genotypetemperature interaction, at 30 °C relative to 20 °C. These 16,916 features were filtered to exclude those without at least 1 Read Per Million (RPM) of 22G-RNA expression in at least three different replicates, resulting in a final set of 14,283 features (13,668 protein-coding genes, 388 pseudogenes, 187 repeats, and 40 transposons) used for further analysis. To assess the similarity of 22G-RNA expression across the 17 samples, we created a multidimensional scaling plot of replicates based on the 22G-RNA expression of these 14,283 genomic features using the plotMDS function from the limma R package (Ritchie et al. 2015). 22G-RNA expression values were normalized to units of log 2 RPM using voom (Law et al. 2014), implemented in the limma R package (Ritchie et al. 2015). We did not include trimmed mean of M-values (TMM) normalization in our normalization pipeline, as our 22G-RNA expression data violate the assumption of TMM that most genes are not differentially expressed, which can lead to inflated false positive rates (Evans et al. 2018). Features with differential 22G-RNA expression between genotype-temperature combinations were detected using limma with the following linear model: 22G-RNA expression = genotype effect + temperature effect + genotype:temperature interaction effect. Features with an adjusted P-value below a falsediscovery rate of 0.05 were considered significant. We then classified features into five non-overlapping groups based on the presence of significant genotype, temperature, or interaction effects: 1) 'genotype only' features have a significant genotype effect, but no significant temperature effect or interaction, 2) 'temperature only' features have a significant temperature effect, but no significant genotype effect or interaction, 3) 'additive genotype-temperature' features have significant genotype and temperature effects, but no significant interaction, 4) 'genotype-temperature interaction' features have a significant interaction effect, and 5) 'no effect' features do not have a significant genotype, temperature, or interaction effect. For each of these five sets of features, we plotted voom-normalized small RNA expression across all 17 samples using the ComplexHeatmap R package (Gu et al. 2016). Expression values for each heatmap row were rownormalized using R's scale function (R Core Team 2021) and ordered by hierarchical clustering. The distribution of protein-coding genotype-temperature interaction genes across each chromosome was visualized using the density calculated by R's hist function (R Core Team 2021), using a bin size of 250 kb.

Comparisons With mRNA Expression
We processed mRNA expression data from Supplementary File S1 of Mark et al. (2019) that correspond to each of the 17 C. briggsae replicate samples used here. This dataset contains voom-normalized mRNA expression data for 16,199 different protein-coding genes, along with classifications describing whether each gene shows significant genotype, temperature, or interaction effects on mRNA expression. We compared expression levels for the 12,623 genes with both 22G-RNA and mRNA expression values available for each replicate in all six genotype-temperature combinations (regardless of expression level), using the mean expression of all replicates in a given treatment. We did not further normalize these gene-level values of 22G-RNA and mRNA expression to log 2 Reads Per Kilobase Million (RPKM), as doing so resulted in significant negative correlations between 22G-RNA/mRNA expression and gene length, although our conclusions with respect to mRNA and 22G-RNA expression do not MBE qualitatively change with normalization to log 2 RPKM rather than log 2 RPM (data not shown). For scatterplots of 22G-RNA versus mRNA expression for these genes, we fit a second-order polynomial curve to each plot and first calculated the correlation between 22G-RNA and mRNA expression of all 12,623 of these genes using Spearman's rank correlation coefficient. We then computed the correlation again using only the genes in the top 10% of 22G-RNA expression in each genotype-temperature combination. To test whether these correlation coefficients were significantly different from 0, we used the cor.test function in R (R Core Team 2021). The significance of differences between correlation coefficients was calculated using a z-test after Fisher's z-transformation (Myers and Sirois 2006). The significance of the overlap between the genes showing a genotype-temperature interaction for both 22G-RNA and mRNA expression was calculated using a hypergeometric test and visualized using the VennDiagram R package (Chen and Boutros 2011).

Expression of Small RNA Pathway Proteins and Targets
To associate trends in 22G-RNA expression with the mRNA expression of proteins involved the 22G-RNA pathway, we extracted Argonaute and RdRP mRNA expression profiles for C. briggsae from the Mark et al. (2019) dataset. To assess the significance of the genotype-temperature interactions on mRNA expression for these genes between 14 °C and 20 °C, and between 20 °C and 30 °C, we modeled voomnormalized mRNA expression for each gene using the linear model: mRNA expression = genotype effect + temperature effect + genotype:temperature interaction effect.
Next, we compiled the identity of genes targeted by these proteins from prior studies of small RNA sequencing that used either immunoprecipitation or mutation of the protein of interest. CSR-1 targets, as determined by CSR-1 immunoprecipitation in C. briggsae, were taken from supplementary material table S2 of Tu et al. (2015). For other 22G-RNA pathway proteins (CSR-1 isoforms CSR-1a and CSR-1b, WAGO-1, RRF-1, EGO-1, and HRDE-1), target genes were defined as the C. briggsae orthologs of target genes from experiments conducted in C. elegans, using the C. briggsae-C. elegans orthologs defined in Tu et al. (2015). Given that 79.6% (3851/4839) of CSR-1 targets in C. briggsae are orthologs of CSR-1 targets in C. elegans (Tu et al. 2015), we expect this orthology approach to be valid presuming that other small RNA pathways show similarly high levels of target conservation between species. This orthology-based approach included genes targeted in C. elegans by the CSR-1 isoforms CSR-1a and CSR-1b [taken from 'IP Enrichment (22G reads)' in supplementary material  Buckley et al. (2012) with an RPKM of at least 50]. Genes that show non-dynamic gene expression across development in C. elegans were defined as the genes classified as having 'flat' gene expression in supplementary material table S1 of Hendriks et al. (2014), again using the C. briggsae orthologs of these C. elegans genes. Wilcoxon rank sum tests were used to calculate the significance of differences in voom-normalized 22G-RNA and mRNA expression of target genes between temperature treatments of the same genotype.

Supplementary Material
Supplementary data are available at Molecular Biology and Evolution online.