Evidence for miRNA‐mediated modulation of the host transcriptome in cnidarian–dinoflagellate symbiosis

Reef‐building corals and other cnidarians living in symbiotic relationships with intracellular, photosynthetic dinoflagellates in the genus Symbiodinium undergo transcriptomic changes during infection with the algae and maintenance of the endosymbiont population. However, the precise regulatory mechanisms modulating the host transcriptome are unknown. Here, we report apparent post‐transcriptional gene regulation by miRNAs in the sea anemone Aiptasia, a model system for cnidarian–dinoflagellate endosymbiosis. Aiptasia encodes mainly species‐specific miRNAs, and there appears to have been recent differentiation within the Aiptasia genome of miRNAs that are commonly conserved among anthozoan cnidarians. Analysis of miRNA expression showed that both conserved and species‐specific miRNAs are differentially expressed in response to endosymbiont infection. Using cross‐linking immunoprecipitation of Argonaute, the central protein of the miRNA‐induced silencing complex, we identified miRNA binding sites on a transcriptome‐wide scale and found that the targets of the miRNAs regulated in response to symbiosis include genes previously implicated in biological processes related to Symbiodinium infection. Our study shows that cnidarian miRNAs recognize their mRNA targets via high‐complementarity target binding and suggests that miRNA‐mediated modulations of genes and pathways are important during the onset and maintenance of cnidarian–dinoflagellate endosymbiosis.


| INTRODUCTION
The trophic and structural basis of one of the most important marine ecosystems, coral reefs, relies on a functional endosymbiosis between cnidarian coral hosts and photosynthetic dinoflagellate symbionts in the genus Symbiodinium, which reside within vesicles (symbiosomes) in the host's gastrodermal cells (Davy, Allemand, & Weis, 2012). In this mutualistic relationship, the host offers a sheltered environment to the alga and provides the inorganic nutrients needed for photosynthesis and growth, whereas the symbiont transfers the majority of its photosynthetic products to the host (Muscatine & Porter, 1977). This tight nutrient coupling allows the animal-alga pair to grow in extremely oligotrophic waters. Most symbiotic cnidarians acquire their endosymbionts with high specificity from the environment ("horizontally") (Dubinsky & Stambler, 2010), either during the larval phase after sexual reproduction or during the recovery from a stress-induced symbiosis breakdown ("bleaching"). There is growing evidence that the establishment and maintenance of the endosymbiosis are accompanied by specific changes in host gene expression, including genes involved in pattern recognition and innate immunity, endocytosis and transmembrane nutrient transport (Baumgarten et al., 2015;Ganot et al., 2011;Lehnert et al., 2014;Rodriguez-Lanetty, Phillips, & Weis, 2006;Schnitzler & Weis, 2010;Voolstra et al., 2009). However, little is yet known about the higher level regulatory mechanisms that might participate in the modulation and orchestration of such changes in gene expression.
MicroRNAs (miRNAs) are small (~22 nt) noncoding RNAs that post-transcriptionally downregulate gene expression through complementary binding of their target mRNAs. In this process, the mature miRNA first binds to the Argonaute (Ago) protein and guides the miRNA-induced silencing complex (miRISC) to partially complementary target sequences within mRNAs, resulting in translational inhibition or target degradation (Carthew & Sontheimer, 2009). Even though the general mechanism of RNA interference and miRNA action is evolutionarily conserved throughout large parts of the phylogenetic tree, several aspects of miRNA biogenesis and mode of action are distinct, in particular between animals and plants (Shabalina & Koonin, 2008). For instance, plant miRNAs bind their mRNA targets mostly with high complementarity within the coding sequence (CDS), leading to target cleavage. In contrast, most metazoan miRNAs only require short complementarity of 2-8 nt in the miRNA 5 0 region (i.e., the "seed") for target recognition in the 3 0 UTR that commonly leads to translational repression (Shabalina & Koonin, 2008). Interestingly, dissection of miRNA function in the "basal" metazoan Nematostella vectensis showed that although miRNA expression led to translational repression (Mauri et al., 2017), the miRNAs also seem to bind their mRNA targets with high, plantlike complementarity leading to mRNA cleavage (Moran et al., 2014).
One miRNA can have a multitude of mRNA targets and widespread downstream effects, making miRNAs regulatory hubs for the fine-tuning of gene expression and the orchestration of broad transcriptomic changes (Esau et al., 2006). Recent studies have suggested that miRNAs are involved in the interactions of hosts with eukaryotic endoparasites such as Leishmania donovani and Plasmodium falciparum (Ghosh, Bose, Roy, & Bhattacharyya, 2013;LaMonte et al., 2012), suggesting that miRNA expression and regulation of target genes might also be involved in transcriptome modulation in the cnidarian-dinoflagellate endosymbiotic relationship.
Using the small sea anemone Aiptasia as a model system for cnidarian-dinoflagellate symbiosis, we took advantage of its recently sequenced genome (Baumgarten et al., 2015) to investigate the potential role of host miRNA regulation by integrating small RNA and mRNA sequencing at different stages of endosymbiont infection.
By establishing a protocol for high-throughput sequencing-crosslinking immunoprecipitation (HITS-CLIP) of the Aiptasia Ago protein, we identified high-confidence ternary interactions involving the Ago effector protein, the potentially regulatory miRNAs, and their target mRNA transcripts on a transcriptome-wide scale. Data produced following this approach supported the conclusion that cnidarian miR-NAs recognize their mRNA targets via high-complementarity target binding and revealed possible miRNA-mediated modulations of genes and pathways important for a functional symbiosis.

| Aiptasia culture
Anemones of the clonal Aiptasia strain CC7 (Sunagawa et al., 2009) were used for all experiments and were cultured as described previously (Baumgarten et al., 2015). Briefly, anemones were kept in autoclaved sea water (ASW) at 25°C on a 12-hr:12-hr light-dark cycle (20-40 lmol photons m À2 s À1 during light cycle) and fed once a week with freshly hatched Artemia (brine shrimp) nauplii larvae.
Water changes were performed on the days following feedings.
Aposymbiotic (dinoflagellate-free) anemones were generated as described previously (Baumgarten et al., 2015) using a combination of cold shocks in the dark at 4°C, treatment with the photosynthesis inhibitor diuron (DCMU) in the light at 25°C, and incubation for ≥6 months in the dark at 25°C with feeding and water changes as before. Prior to experimental use, aposymbiotic Aiptasia were acclimatized to a 12-hr:12-hr light-dark cycle for at least two weeks at 25°C and checked by fluorescence stereomicroscopy to verify that they were still devoid of dinoflagellates (as judged by the absence of chlorophyll autofluorescence).

| Infection with Symbiodinium and RNA sequencing
To analyse RNA from anemones at different stages of infection, aposymbiotic anemones were infected with the clonal, axenic Symbiodinium minutum strain SSB01 (Baumgarten et al., 2015;Xiang, Hambleton, DeNofrio, Pringle, & Grossman, 2013) essentially as described previously (Baumgarten et al., 2015). Briefly, anemones of 5-10 mm oral-disc diameter were held in 6-well plates (Falcon # 353046); each well contained 3-5 anemones and 15 ml of autoclaved and sterile-filtered sea water (ASFW). After 16 d of acclimatization (see above), infection proceeded as follows: day 1, addition of Symbiodinium cells (from a culture growing exponentially in F/2-Si medium; Baumgarten et al., 2013) to a final concentration of~10 5 cells/ml; day 2, feeding (as above) without change of ASFW or adding additional algae; day 3, change of ASFW and addition of algae at~10 5 cells/ml; day 11, end of incubation with Symbiodinium cells by moving the anemones into ASFW without algae in fresh 6-well plates; and days 11-30, regular feeding and ASFW changes as above. Samples were taken for RNA analysis at days 1 (before addition of algae; aposymbiotic anemones), 12 (partially populated anemones), and 30 (by which time the anemones should be fully populated) (Hambleton, Guse, & Pringle, 2014 Table S1). The processing of mRNA samples and the analysis of differentially expressed genes from the infection experiment were described by (Baumgarten et al., 2015); the mRNA sequences can be accessed at the NCBI SRA (Accession nos.: aposymbiotic, SRX757525; partially populated, SRX757526; symbiotic, SRX757528). To estimate the progression of Symbiodinium infection over the 30-day time-course, the mRNA sequences were mapped to a S. minutum reference transcriptome (NCBI BioProject Accession no. PRJNA274852) using bowtie2 (Langmead & Salzberg, 2012). The percentages of read pairs that aligned to the symbiont transcriptome without mismatches are shown in Figure S1a.

| miRNA annotation
The small RNA reads were first trimmed by removing adapter sequences and filtered to remove low-quality reads using Trimmomatic (Bolger, Lohse, & Usadel, 2014). Only reads of 15-32 nucleotides were retained, and redundant sequences were collapsed using the collapse_reads.pl script from the miRDeep2 package (Friedlander, Mackowiak, Li, Chen, & Rajewsky, 2012). Next, we annotated the regions of the Aiptasia genome encoding noncoding RNAs (ncRNAs), including snoRNAs, snRNAs, tRNAs and rRNAs using reference ncRNA sequences retrieved from the RNAcentral database (The RNAcentral Consortium, 2015) and aligned to the Aiptasia genome using BLASTN (E-value cut-off ≤1e-5) (Altschul, Gish, Miller, Myers, & Lipman, 1990). The Aiptasia ncRNA regions identified by this alignment were then used as a database for mapping of the small RNA reads using bowtie2 (Langmead & Salzberg, 2012). Small RNA reads that mapped to the ncRNA database using bowtie2 were removed before further analyses. miRNAs were then annotated using the miRDeep2 package with default settings (Friedlander et al., 2012). To identify putatively conserved miRNAs based on previous de novo annotations of other cnidarian genomes, we created a reference library of mature miRNA sequences from N. vectensis (Grimson et al., 2008;Moran et al., 2014), Hydra magnipapillata (Krishna et al., 2013) and Stylophora pistillata (Liew et al., 2014). An initial set of 145 miRNAs was predicted and further inspected for characteristics of bona fide miRNAs. These included (i) potential folding of the precursor miRNA transcript (pre-miRNA) into an unbranched hairpin structure, (ii) 5 0 ends of small RNA reads mapping uniformly to the pre-miRNA hairpin, (iii) differing frequencies of mature miRNA and star miRNA (i.e., the complementary strand that is not bound to the Ago protein) reads mapping to the pre-miRNA stem and (iv) 2-nt overhangs of the mature miRNA and the star miRNA sequence in the miRNA duplex (sensu Baumgarten et al., 2013;Tarver, Donoghue, & Peterson, 2012). Applying these additional stringent criteria, we identified 46 miRNAs mapping to 60 genomic loci (Table S2). We then used the built-in scoring mechanism of miRDeep2 to assign probabilities of the 46 miRNA to be true positives, performing 100 rounds of permutation (Friedlander et al., 2008(Friedlander et al., , 2012. All 46 miRNAs had true-positive probabilities of detection ≥84%, except for the evolutionarily conserved apa-mir-9454. Reviewing the prediction for apa-mir-9454 showed that no star sequences for this miRNA were recovered in the Aiptasia small RNA libraries, decreasing its probability of identification by the algorithm. However, because of the perfect conservation of the mature miRNA sequence, we retained apa-mir-9454 for subsequent analyses. 2.4 | miRNA differential expression and motif analysis To assess expression levels of miRNAs, we mapped each set of preprocessed small RNA reads (i.e., small RNAs retained after ncRNA filtering) independently to the identified miRNA precursor sequences using the script quantifier.pl from the miRDeep2 package (Friedlander et al., 2012). The miRNA read counts of each replicate were then normalized to the total number of mapped small RNA reads to allow for comparison between samples. These analyses were performed with three replicates for each time point of the infection experiment (see above) using the EDGER package (Robinson, McCarthy, & Smyth, 2010)   To identify motifs putatively involved in the transcriptional control of the one apa-mir-2022a and two apa-mir-2022b loci (see Section 3), we extracted the 1 kb of sequence upstream of the pre-miRNA start position of each locus. The motif analysis was performed using the MEME web suite (Bailey et al., 2009), searching for up to 10 motifs that occur zero or only one time per region. For comparison, the 1 kb of sequence downstream of each locus was also used for motif identification using the same settings.
2.5 | Generation and testing of anti-AipAgo1 antibody Aiptasia Argonaute (AipAgo) proteins were identified by BLASTP alignments of the N. vectensis Ago protein sequences (UniProt Accession nos. U3MIG1, U3MH35) to the Aiptasia genomic protein set (http://aiptasia.reefgenomics.org) (E-value ≤1e-5). From the canonical AipAgo1, a peptide of 25 amino acids (residues 811-835) ( Figure S1b) was selected to be used as antigen. To ensure that it was specific to AipAgo1, the peptide sequence was searched against the entire Aiptasia genomic protein set (Baumgarten et al., 2015) using BLASTP (E-value cut-off ≤1e-5) (Altschul et al., 1990); no significant alignment was found other than the corresponding position in AipAgo1 itself. This peptide was subsequently synthesized and used for the generation of polyclonal antibodies (Acris Antibodies, Herford, Germany). Briefly, antibodies (Ab) were raised in rabbits following immunization with 1 mg of KLH (keyhole limpet haemocyanin)-conjugated peptide antigen. The antigen-specific Ab titre was determined with an enzyme-linked immunosorbent assay (ELISA) using horseradish peroxidase (HRP) conjugated to anti-rabbit-IgG Ab and ATBS (2,2 0 -azinobis [3-ethylbenzothiazoline-6-sulfonic acid]-diammonium salt) as substrate. Abs were affinity-purified from the serum of the highest titre bleed using the synthesized peptide antigen on a cyanogen bromide-activated agarose column.
To test the specificity of the antibodies by Western blotting, Aiptasia lysates were prepared by sampling~15 adult anemones (5-10 mm oral-disc diameter) in ice-cold PBS and grinding them into pieces of several mm 3 with a plastic mortar and pestle. Following centrifugation (200g, 10 min, 4°C), the tissue pellet was resuspended in three volumes of NP-40 cell lysis buffer (Thermo Fisher # FNN0021) supplemented with 1 mM PMSF and a protease-inhibitor cocktail (Sigma-Aldrich # P2714) according to the manufacturers' instructions. The tissue samples were then disrupted using a Wheaton glass homogenizer, incubated on ice for 30 min and centrifuged (16,000g, 10 min, 4°C); the protein-containing supernatant was then sampled for further processing. Proteins were resolved electrophoretically on a 12% SDS-polyacrylamide gel for~1.5 hr at 120 V using a 2-gel Tetra and Blotting Module (Bio-Rad # 1660828EDU) and transferred to a PVDF membrane for 2 hr at 100 V using the same module. After blocking overnight at 4°C using StartingBlock T20 (PBS) (Thermo Fisher # 37538), the membrane was incubated for 2 hr at 37°C with 20 lg of Ab diluted 1:1,000 in Tris-buffered saline (pH = 7.5) supplemented with 0.05% (v/v) Tween-20 (TBST), followed by 3-hr incubation at room temperature with a goat-anti-rabbit-IgG secondary Ab tagged with DyLight Fluor 488 (Thermo Fisher # A-11034) diluted 1:10,000 in TBST. Each incubation with antibodies was followed by extensive washes with TBST to remove any nonspecific binding. Protein bands were visualized using a Typhoon Scanner 9410 (GE Healthcare).

| Immunohistochemical staining
For whole-mount staining, anemones were placed in a six-well plate with ASFW, stunned for 15 min with 0.2% magnesium chloride and then incubated in 1% formaldehyde for 20 min at room temperature.
The fixed anemones were transferred to 1.5-ml microcentrifuge tubes and washed in 500 ll PBST (phosphate-buffered saline with 0.1% Triton X-100) for 15 min on a spin-wheel at room temperature.
For antigen exposure, the PBST was removed, and the anemones were suspended in 500 ll of 10 nM sodium citrate in sterile water and incubated in a water bath at 99°C for 30 min. Subsequently, the anemones were washed with PBST and incubated for 2 hr in PBSTB blocking solution at room temperature. For staining, the samples were incubated in anti-AipAgo1 primary antibody (dilution: 1:250 in PBSTB) or in PBSTB alone (negative control) overnight at 4°C on a spin-wheel. The samples were washed twice with PBSTB for 10 min per wash, followed by incubation with Dylight Fluor 650-tagged goat-anti-rabbit-IgG secondary antibody (Thermo Fisher #84546) for 2 hr at room temperature. DNA staining followed the protocol described below.
Aposymbiotic anemones were macerated according to (David, 1973). Briefly, whole animals were placed in a maceration solution containing glycerine, glacial acetic acid and water (ratio 1:1:13) at room temperature. The tissue was left soaking for 30 min with F I G U R E 1 Overview of the AipAgo1 cross-linking immunoprecipitation (CLIP) experiment and sequence analysis. (a) Workflow for the CLIP experiments (see Section 2). The tripartite miRNA-mRNA-Ago complex was cross-linked with UV light in vivo followed by cell lysis and immunoprecipitation of AipAgo1 and the bound RNA. After digestion of mRNA with RNase A/T1 (preserving only tags protected by AipAgo1 around the miRNA binding site), the protein-RNA complexes were resolved by molecular size using SDS-gel electrophoresis and blotted to a nitrocellulose membrane. The region corresponding to the molecular weight of AipAgo1 (~96 kDa) was excised (Panel c), protein was digested with proteinase K, and RNA was extracted and prepared for sequencing (Panel e). TSAP, thermosensitive alkaline phosphatase; T4 PNK, T4 polynucleotide kinase. (b) Red line, electropherogram showing length distribution of RNA tags extracted immediately after step 3 (Panel a; see Section 2). Blue line, same with a control IgG; no RNA tags were detected after immunoprecipitation and RNase digestion. (c) Western blot of AipAgo1-RNA complexes. The red box indicates the membrane area that was excised for further processing.    Beads diluted in the premixed elution solution were heated for 10 min at 70°C, and the supernatant was collected for SDS-PAGE analysis (see below).
As an initial positive control for successful RNA co-immunoprecipitation, RNAs were isolated directly from the Dynabeads (after RNase digestion and washing) using TRIzol LS (Thermo Fisher # 10296010) and the mirVana small-RNA-enrichment kit (Thermo Fisher # AM1560) and visualized on a Bioanalyzer 2100 (Agilent Technologies, RNA Pico Chip) ( Figure 1b). As a negative control, a mock immunoprecipitation was performed using a nonspecific goat anti-mouse-IgG Ab (Thermo Fisher # 35502), and RNAs that might have copurified with the antibody were again isolated directly from the beads and analysed as just described (Figure 1b). RNAs were then isolated from the regions of the nitrocellulose membranes that should contain the AipAgo1 protein (predicted molecular weight~96 kDa) and larger RNA -AipAgo1 complexes ( Figure 1c). The membrane regions were excised, cut into 1-to 2mm squares and collected in RNase-free microcentrifuge tubes.

| HITS-CLIP sequence analyses
As outlined in Figure 1e, sequences obtained from the HITS-CLIP experiment were first trimmed of adapter sequences and filtered to remove low-quality sequences using Trimmomatic (Bolger et al., 2014). For the detection of miRNAs, reads were then aligned to the pre-miRNA sequences obtained from the miRNA annotation (miRDeep2) using NovoAlign (http://www.novocraft.com), requiring a minimum alignment of 15 nucleotides and allowing soft clipping of the read ends (Figure 1e, left). The numbers of reads mapping to each pre-miRNA were then tested for correlation between the replicates for each biological sample after first normalizing the raw counts for each sample by the total number of CLIP sequences that aligned to pre-miRNAs in the respective library (Figure 1e (Table S3). Overlapping CLIP-tags (i.e., CLIP-tag peaks) and their respective heights (i.e., numbers of CLIP-tags per peak) were then identified, requiring a minimum overlap of one nucleotide within a peak. The bedtools "intersect" tool (Quinlan, 2014) was used to identify genomic regions that showed CLIP-tag peaks across all samples.
The CLIP-tag heights were normalized as for the CLIP miRNA counts (see above), and the normalized heights of the CLIP-tag peaks containing three or more CLIP-tags at the same genomic region in each sample were tested for correlation.
Pyicoclip ( were then extracted from the Aiptasia genome using bedtools "getfasta" (Quinlan, 2014). Next, RNAhybrid (Kr€ uger & Rehmsmeier, 2006) was used to map the 46 bona fide miRNAs (see above) to the CLIP-tag peak sequences requiring (i) no a priori helix constraints and (ii) a minimum free energy of binding (MFE) ≤À15 kcal/mol, while allowing (iii) GU wobbles in the miRNA-target alignment and (iv) a maximum of 1-nt bulges on either side of the alignment. When a CLIP-tag peak region was targeted by multiple miRNAs, we only retained the interaction that showed the highest number of binding nucleotides between the miRNA and the CLIP-tag peak region.
For comparison of miRNA-mRNA complementarity between Aiptasia and published mammalian CLIP data sets, we retrieved the human miRNA CLIP data that were initially sampled by Kishore et al. (2011) and reanalysed by Grosswendt et al. (2014) for the presence of chimeric miRNA-mRNA reads to unambiguously identify miRNA-mRNA target sites. We mapped the human miRNA to their respective mRNA target site using RNAhybrid (Kr€ uger & Rehmsmeier, 2006) and calculated the numbers of interacting miRNA nucleotides with the same settings as described above.

| Annotation of conserved and species-specific miRNAs in Aiptasia
We investigated the Aiptasia miRNA repertoire using replicate small RNA libraries prepared from Aiptasia polyps at different stages of infection (aposymbiotic, partially populated and fully populated) with the endosymbiotic S. minutum strain SSB01 (see Section 2, Figure 2a and Figure S1a, and Table S1). Initial inspection of the small RNA reads revealed a length distribution with an apparent overrepresen- From the pooled small RNA samples, we annotated 46 high-confidence miRNAs, transcribed from a total of 60 distinct genomic loci (Table S2) . At each length, the small RNA reads are sorted by their 5 0 nucleotide (coloured bars), and the sum at each length is shown as the grey area. (c) Example of a miRNA (apa-mir-2022a) matching all criteria for bona fide miRNA annotation. Box, predicted secondary hairpin structure of the pre-miRNA as predicted using RNAfold (Hofacker, 2003) within miRDeep2 (Friedlander et al., 2012). Bottom, the relative frequencies of small RNA reads that map to the opposite sides of the hairpin (ordinate) are shown along an axis corresponding to the primary pre-miRNA sequence (abscissa). Reads matching the mature miRNA sequence (red) have homogeneous 5 0 ends and considerably higher read counts than do those matching the complementary strand (star-miRNA; purple). Both star and mature miRNA sequence have a distinct length of 23 nt. (d) Overlap of evolutionarily conserved mature miRNAs in the anthozoans for which data are available; the numbers of conserved and distinct species are indicated. Only six miRNAs (denoted by the *) are found in all four species: mir-100, mir-2022, mir-2023, mir-2030, mir-2036 and mir-2050 (Table S3). Examination of the clusters of overlapping CLIPtags (≥1-nt overlap) at distinct genomic loci showed that the numbers of CLIP-tags per cluster also correlated well both between the two replicates for each symbiotic state and between the two symbiotic states (Figure 3c).
The high correlation between replicates observed with both the Ago-miRNA (Figure 3d) and Ago-mRNA (Figure 3e)   . Note accumulation of AipAgo1 in particular locations that presumably correspond to P-bodies. (d) Correlations between the miRNA read counts (log10 values) from the independent cross-linking immunoprecipitation (CLIP) experiments for the two replicates of aposymbiotic anemones (left), the two replicates of fully populated symbiotic anemones (middle), and the means of the two replicates for each of the two symbiotic states (right). (e) Correlation of the CLIP-tag peak heights across the genome in the two independent replicates of the CLIP experiment with aposymbiotic animals (A 1 , A 2 ) and the two replicates with symbiotic animals (S 1 , S 2 ). The rows represent the various genomic regions with ≥3 CLIP-tags in each of the replicates; the columns represent the four replicates. The correlation coefficients are shown below. (f) Frequencies (%) of miRNA-mRNA interactions with various numbers of miRNA nucleotides involved in mRNA target binding (blue area: Aiptasia CLIP data [this study]; grey area: Human CLIP data; Grosswendt et al., 2014;Kishore et al., 2011, left-hand axis). The average minimum free energy of binding for each of the miRNA-mRNA interactions in the Aiptasia data is shown as orange diamonds (right-hand axis). (g) Genome-browser view of the region containing a gene (systematic name AIPGENE22434) that encodes a cnidarian Ficolin-like protein (CniFL; see Baumgarten et al., 2015). Most CLIP-tags of the pooled aposymbiotic and symbiotic samples mapped to a single narrow region (blue boxes below Exon 1 in the gene model), which putatively represents a target site of apa-mir-12452. A few other CLIP-tags mapped to the 3 0 end of Exon 3; the significance of this mapping is unclear.
The box at lower left shows the sequence of apa-mir-12452 complementary to the putative mRNA target site. Green, untranslated regions; yellow, coding sequences; narrow lines, introns; arrows, directions of transcription [Colour figure can be viewed at wileyonlinelibrary.com] robustness, as observed also in previous studies (Chi et al., 2009;Weyn-Vanhentenryck et al., 2014). In addition, the interactions identified in the aposymbiotic and symbiotic samples (Figure 3d,e) revealed only limited differences in the overall numbers of interactions. The fine-tuned regulation of mRNA transcripts levels might therefore originate from quantitative (i.e., differential expression of targeting miRNAs) rather than from qualitative (i.e., different miRNA-mRNA interactions) differences. Thus, although it is possible that infection state-specific differences in miRNA-mRNA-Ago inter- actions may yet be detected, we pooled the CLIP-tags of all four biological samples for subsequent analyses.
We next associated significant CLIP-tag peaks with exonic regions using the modified false discovery rate (FDR) algorithm proposed by Yeo et al. (2009)  of these peaks lay within the protein-coding sequences of transcripts, with an additional 27% and 10% falling within 5 0 and 3 0 UTRs, respectively.
To identify functionally significant sites of miRNA-mRNA binding, we next mapped the bona fide mature miRNAs annotated by miRDeep2 to the identified CLIP-tag peak sequences (i.e., the genomic sequence from the start to the end of the CLIP-tag peak). This analysis identified 3,377 interactions between a miRNA and an mRNA-Ago site with a minimum free energy of binding of ≤À15 kcal/mol (Table S4). A total of 2,758 of these interactions included a single CLIP-tag peak sequence and a single miRNA, whereas 619 CLIP-tag peak sequences showed more than one equally good mRNA interaction with a miRNA. Further inspection of these 619 interactions showed that 226 CLIP-tag peak sequences were targeted by two almost identical miRNAs, apamir-2022a and apa-mir-2022b, and another 30 CLIP-tag peak sequences were targeted by a single miRNA multiple times and equally well (i.e., with the same number of binding miRNA nucleotides; see Section 2).

(a) (c) (b)
F I G U R E 4 Differential expression of miRNAs during establishment of the Symbiodinium endosymbiosis. (a) Heatmap of miRNAs with significant differential expression (FDR ≤ 0.01); most are upregulated in the symbiotic relative to the aposymbiotic state. The columns represent the six independent experiments (three with aposymbiotic animals and three with fully populated symbiotic animals). (b) Top, relative positions of the apa-mir-2022b and apa-mir-2022a loci in the Aiptasia genome. Grey arrows, distances in kilobases (kb) between the miRNA loci; red arrows, directions of transcription. Bottom, presence of apa-mir-2022a, but not apa-mir-2022b, also in other anthozoan genomes, the single nucleotide difference between the two miRNAs is highlighted in red. (c) Locations of conserved sequence motifs detected in the regions 1-kb upstream of the apa-mir-2022a and apa-mir-2022b loci (top) but not in the regions 1-kb downstream of the loci (bottom); in the downstream regions, only short fragments of the sequence motifs are present, and in seemingly random order. Each coloured box denotes a sequence motif that is conserved between two or among all three of the genomic loci (see Figure S3 for details of the actual sequences) [Colour figure can be viewed at wileyonlinelibrary.com] In mammals, it has been shown that binding of as few as seven nucleotides in the 5 0 seed region of a mature miRNA is sufficient for effective miRNA-mRNA binding (Bartel, 2009). In contrast, many of the interactions identified here seemed to include a much higher number of binding nucleotides from the miRNA, especially when directly compared to human CLIP data (Grosswendt et al., 2014;Kishore et al., 2011) (Figure 3f). While most human miRNA-mRNA interactions involve only eight miRNA nucleotides, nearly a third of the Aiptasia miRNA-mRNA interactions (n = 1,061) have a complementarity to the mRNA target of ≥17 miRNA nucleotides (Figure 3f).
Such high complementarity has also been observed in other cnidarians (Moran et al., 2014)  To explore the possible coregulation of these three miRNAs, we searched for conserved sequence motifs in the 1-kb regions upstream of their loci. Interestingly, we found that the two apa-mir-2022b loci, separated by~275 kb, share several long and highly conserved motifs in their upstream regions, whereas the upstream region of apa-mir-2022a is distinct (Figure 4c, top and Figure S3).
The regions downstream of the three loci also show no comparable sequence conservation (Figure 4c, bottom). Because the exact transcription start sites of the primary (pri-) miRNA transcripts for the three loci are unknown, the motifs conserved between the two apamir-2022b loci could also lie within the pri-miRNA transcripts, so that the sequence conservation observed could be involved in cocontrol of either transcription or post-transcriptional processing of the pri-miRNA transcripts.

| Possible modulation of genes involved in symbiont acquisition and/or maintenance by Aiptasia miRNAs
Post-transcriptional control of gene expression by miRNAs usually results in the downregulation of cognate target genes either through endonucleolytic cleavage or inhibition of mRNA translation (Ghildiyal & Zamore, 2009). To assess possible changes of mRNA transcript levels due to the interactions with miRNAs, we measured mRNA levels in the same total RNA samples from the infection experiment that were used for the miRNA-expression and CLIP-tag analyses. In total, we found 885 interactions identified by CLIP between mRNA transcripts and the 12 miRNAs that showed significant differential expression, and thus are more likely to have a regulatory role in the endosymbiosis (see above).
For the 11 differentially expressed miRNAs that were upregulated in the symbiotic state (Figure 4a), we found interactions with a total of 289 putative target genes that also showed higher expression in the symbiotic state (log2 fold change > 0), 228 putative target genes that showed lower expression in the symbiotic state (log2 fold change < 0), and 108 putative target genes with no change in transcript abundance. For the single differentially expressed miRNA that was downregulated in the symbiotic state (apa-mir-2026), we found seven putative target genes that were upregulated and four that were downregulated, along with five putative target genes that did not show a measurable change in transcript abundance (Table S4). Thus, there was no general correlation in which changes in miRNA levels led consistently to inverse changes in the levels of the predicted target mRNAs.
However, upon closer examination of the putative target genes, we found several whose products have previously been implicated in the onset and/or maintenance of cnidarian-dinoflagellate endosymbiosis. These include a membrane receptor of the transforming growth factor b receptor family (TGFbR) (Detournay, Schnitzler, Poole, & Weis, 2012), in whose mRNA three sites were targeted by apa-mir-12450, apa-mir-12428, and apa-mir-2037, and the downstream messenger of TNF receptors, TRAF (Barshis et al., 2013), targeted at one site in its mRNA by apa-mir-12435. We also identified components of a putative fibroblast-growth-factor (FGF) signalling cascade including the extracellular ligand FGF (two sites, targeted by apa-mir-12424 and apa-mir-12448), its receptor (FGFR, five sites targeted by apa-mir-12449, apa-mir-12421, apa-mir-12438, apa-mir-BAUMGARTEN ET AL. | 413 12441, apa-mir-2022a and apa-mir-2025), and a downstream intracellular messenger, GRB2 (one site, targeted by apa-mir-2022a and apa-mir-2022b), to be under possible post-transcriptional miRNA control. Although we saw no changes in levels of the TGFbR or GRB2 mRNAs, we found that the FGF, TRAF and FGFR mRNA decreased considerably in the symbiotic state, possibly reflecting the increases in the corresponding miRNAs (Table S4).

| DISCUSSION
The establishment and maintenance of cnidarian-dinoflagellate endosymbiosis is accompanied by marked changes in host gene expression (Baumgarten et al., 2015;Ganot et al., 2011;Lehnert et al., 2014;Rodriguez-Lanetty et al., 2006;Schnitzler & Weis, 2010;Voolstra et al., 2009). In this study, we identified the miRNAs of Aiptasia and analysed their expression patterns and mRNA targets in order to illuminate their possible roles in regulating these changes in gene expression. To do so, we first established high-throughputsequencing cross-linking immunoprecipitation (HITS-CLIP) in Aiptasia.
This protocol can be applied readily to any other cnidarian as well as any RNA-binding protein and its target mRNA transcripts of interest.  (Darnell, 2013;Licatalosi et al., 2008;Ule et al., 2003) and endoparasite pathogenesis (Vembar, Macpherson, Sismeiro, Copp ee, & Scherf, 2015). More broadly, we hope that the protocol described here will provide a foundation for further studies of RNA-protein interactions in basal metazoans in general and in particular in regard to their putative role in cnidarian-dinoflagellate endosymbiosis.
Our annotation of the Aiptasia miRNAs revealed a surprisingly indicate lineage-specific functional divergence. Moreover, the conservation of an individual nucleotide substitution in the 3 0 region of the two apa-mir-2022b copies compared to apa-miR-2022a suggests a functional importance for nucleotide binding beyond the 5 0 "seed" binding. This interpretation is consistent with the generally large hybridization length of the miRNA-mRNA duplex (i.e., the number of miRNA nucleotides involved in mRNA-target binding) observed in our CLIP data, in particular when compared to human miRNA-mRNA interactions. Similarly, miRNA-target studies in N. vectensis showed that besides miRNA-mediated translational repression (Mauri et al., 2017), cnidarian miRNAs might often lead also to target cleavage, presumably due to much higher complementarity between the miRNA and mRNA target transcript, as seen also in plant miRNAs (Moran et al., 2014).
Host gene expression is probably regulated by several mechanisms in response to Symbiodinium infection. In this context, the lack of a clear correlation between miRNA and target mRNA expression might point towards a miRNA-mediated, fine-scale modulation rather than a global regulation of gene expression. Furthermore, certain miRNA-mRNA interactions and the resulting putative modulation of target transcript levels might be restricted to a particular set of cells or tissues and so would not have been resolved in the current study.
Even though the identified interactions cannot currently be functionally validated due to the absence of appropriate techniques to knock down or overexpress miRNAs in Aiptasia, the significant changes in miRNA expression do suggest that miRNA-mediated modulation of host gene expression constitutes one mechanism regulating the endosymbiotic relationship.
In particular, the identification of differentially expressed miRNAs that interact with mRNAs whose products are commonly considered to play key roles in endosymbiosis onset and/or maintenance can provide a model of how miRNA action may be involved in this process ( Figure 5). For instance, the interaction of miRNAs with the mRNAs for FGFR, TGFbR and components of the TNFR/TRAF pathways ( Figure 5) is particularly intriguing, as these proteins were previously suggested to be involved in the priming of immune and stress responses of symbiotic cnidarians, pathways that might have been co-opted during the evolution of mechanisms allowing a stable endosymbiosis (Barshis et al., 2013;Detournay et al., 2012). In addition, proteins that localize to the symbiosome membrane are thought to be critical for the establishment and/or maintenance of the endosymbiosis, including proteins involved in symbiosome maturation (Chen, Cheng, Hong, & Fang, 2004;Chen, Cheng, Sung, Kuo, & Fang, 2003) and in cross-membrane nutrient exchange (Dani, Ganot, Priouzeau, Furla, & Sabourault, 2014;Dani et al., 2017;Lehnert et al., 2014). In this context, we also identified miRNA-interacting mRNAs that encode proteins that are likely to be involved in the maintenance of the symbiosome, including a homolog of LAMP1, a protein related to lysosome-and potentially symbiosome-trafficking and maturation (Mohamed et al., 2016); a homolog of the sterol transporter NPC1 (Dani et al., 2014(Dani et al., , 2017Ganot et al., 2011;Lehnert et al., 2014); and a homolog of the peptide transporter ABCB9 (Davy et al., 2012). Of note, LAMP1 and NPC1 mRNA levels were reduced in symbiotic anemones, perhaps reflecting the binding of their regulatory miRNAs and arguing for a repressive effect in these particular cases (Table S4).
On a broader scale, recent studies on the interaction of plant roots and fungal pathogens have suggested that small RNAs (sRNAs) might also act in cross-kingdom RNA interference and regulation of gene expression (Weiberg, Wang, Bellinger, & Jin, 2014). Here, sRNAs of the fungal pathogen are expressed and translocated to the host cell and bound to plant Argonaute proteins to specifically downregulate host immunity genes (Weiberg et al., 2013). Even though we did not find previously identified Symbiodinium miRNAs or sRNAs that were bound to AipAgo1 in the experiments described here, a similar mechanism was recently hypothesized to function also in the coral-Symbiodinium endosymbiosis (Baumgarten et al., 2013;Lin et al., 2015).
Taken together, our data suggest that miRNAs are involved not only in defining the susceptibility of host organisms to endoparasites (Ghosh et al., 2013;LaMonte et al., 2012;Weiberg et al., 2014) but also in post-transcriptional modulations of gene expression that might be important in mutualistic relationships such as the dinoflagellate endosymbiosis in anemones and corals. With techniques now available such as the HITS-CLIP protocol described here, studies of this and other symbiotic systems must-and can-now also incorporate analyses of post-transcriptional levels of regulation.

ACKNOWLEDG EMENTS
Research reported in this publication was supported by the King Abdullah University of Science and Technology (KAUST) and by Gordon and Betty Moore Foundation Grant 2629.01 (to J.R.P.). We thank the Bioscience Core laboratory and the Imaging and Characterization Core laboratory at KAUST for help with the protein characterization and microscopy.

CONFLI CT OF INTERESTS
The authors declare no conflict of interests.