Viral suppressor of RNA silencing in vascular plants also interferes with the development of the bryophyte Physcomitrella patens.

Plant viruses are important pathogens able to overcome plant defense mechanisms using their viral suppressors of RNA silencing (VSR). Small RNA pathways of bryophytes and vascular plants have significant similarities, but little is known about how viruses interact with mosses. This study elucidated the responses of Physcomitrella patens to two different VSRs. We transformed P. patens plants to express VSR P19 from tomato bushy stunt virus and VSR 2b from cucumber mosaic virus, respectively. RNA sequencing and quantitative PCR were used to detect the effects of VSRs on gene expression. Small RNA (sRNA) sequencing was used to estimate the influences of VSRs on the sRNA pool of P. patens. Expression of either VSR-encoding gene caused developmental disorders in P. patens. The transcripts of four different transcription factors (AP2/erf, EREB-11 and two MYBs) accumulated in the P19 lines. sRNA sequencing revealed that VSR P19 significantly changed the microRNA pool in P. patens. Our results suggest that VSR P19 is functional in P. patens and affects the abundance of specific microRNAs interfering with gene expression. The results open new opportunities for using Physcomitrella as an alternative system to study plant-virus interactions.

INTRODUCTION

Phytopathogens, such as viruses and fungi, suppress RNA silencing in plants and hamper plant defense by altering small RNAs in plants (Qiao et al., 2013; Yin et al., 2019). Plant viruses are important pathogens that cause significant economic losses in crop plants and are able to infect and induce mild to severe symptoms in a wide variety of different plant species (K. B. G. Scholthof et al., 2011). Infected plants show symptoms such as leaf yellowing, leaf and growth distortions including curling of the leaf, stunting of the whole plant and abnormalities in flower and fruit formation (Agrios, 2005).

Most plant viruses have single‐stranded, positive‐sense RNA genomes, and their replication occurs through RNA intermediates (Gergerich & Dolja, 2006). Plants recognize virus‐derived molecules and this recognition induces defense responses (Daròs, 2017). RNA silencing, also called RNA interference (RNAi), is a sequence‐specific post‐transcriptional regulatory system with diverse biological roles, including defense against viral infection. It is conserved among plants, including bryophytes (Bezanilla, Pan, & Quatrano, 2003; Burgyan & Havelda, 2011; Hamilton & Baulcombe, 1999; Khraiwesh et al., 2010).

Antiviral RNAi begins with the recognition of viral RNA by Dicer‐like (DCL) proteins, which process double‐stranded viral RNA into different types of viral small interfering RNA duplexes (vsiRNAs) (Garcia‐Ruiz et al., 2010; Qu, Ye, & Morris, 2008; X. B. Wang et al., 2011). vsiRNAs are incorporated into RNA‐induced silencing complexes (RISCs) where one strand (the passenger strand) is removed and the other strand (the guide strand) is retained to direct the RISC to the viral RNA target (Carbonell et al., 2012; Carthew & Sontheimer, 2009; Zhang et al., 2006). Argonaute (AGO) proteins are catalytically active slicers of RISCs and promote the degradation of the viral RNA target.

To overcome antiviral RNA silencing, viruses have evolved viral suppressors of RNA silencing (VSR). Viral suppressors are highly diverse and operate through different mechanisms, for example, by directly interfering with vsiRNA duplexes and preventing RISC formation, or indirectly by interacting with components of the RISC itself or by modulating the plant immune responses in an early stage of infection by regulating microRNA (miRNA) activity (Csorba, Kontra, & Burgyan, 2015; Pertermann et al., 2018). In vascular plants, most VSR‐mediated inhibition of RNAi occurs through sequestration of small RNA (sRNA) duplexes by binding to double‐stranded RNAs (dsRNAs) or through physical interaction with argonaute1 (AGO1) to prevent siRNA or miRNA loading (Burgyan & Havelda, 2011; S. R. Liu, Zhou, Hu, Wei, & Zhang, 2017; M. B. Wang, Masuta, Smith, & Shimura, 2012). The VSR P19 of tombusvirus binds and sequesters many plant sRNAs and vsiRNA to suppress their interaction with AGO1 (Csorba et al., 2015; S. R. Liu et al., 2017). VSR can also interfere with antiviral silencing at multiple points (Csorba et al., 2015). For example, the VSR 2b of cucumber mosaic virus (CMV) has many targets within the antiviral silencing pathway, and it interacts with the PAZ domain of AGO1 to interfere with RISC activity. 2b can also prevent the spread of the long‐range silencing signal and interferes with DNA methylation in the nucleus (Duan et al., 2012; Guo & Ding, 2002). Furthermore, P19 can modulate endogenous miR168 levels to reduce cellular levels of AGO1 in Nicotiana benthamiana plants, indicating that one P19 mechanism alleviates an anti‐viral function of AGO1 (Várallyay, Válóczi, Agyi, Burgyán, & Havelda, 2017).

Considering that viruses are the most abundant biological entities on Earth (Edwards & Rohwer, 2005) and infect all types of organisms (Retel, Märkle, Becks, & Feulner, 2019), surprisingly few studies describe virus interaction or occurrence in bryophytes. Maumus, Epert, Nogue, and Blanc (2014) and Lang et al. (2018) found nucleocytoplasmic large DNA virus (NCLDV) insertions within the genome of P. patens. Moreover, Maumus et al. (2014) found that of 13 studied land plant species (Arabidopsis thaliana, Glycine max, Medicago truncatula, Populus trichocarpa, Ricinus communis, Vitis vinifera, Solanum lycopersicum, Oryza sativa, Sorghum bicolor, Amborella trichopoda, Picea abies, P. patens and Selaginella moellendorffii), two (P. patens and S. moellendorffii) contain NCLDV‐like genes, suggesting that land plants have been infected by NCLDV relatives. Stough et al. (2018) used metatranscriptomes to study the diversity and activity of virus populations in the Sphagnum microbiome. They screened for the presence of conserved virus marker genes and identified 33 contigs representing phages, 114 contigs originating from undescribed single‐stranded RNA viruses and 64 contigs representing NCLDVs. Their results showed that the Sphagnum phyllosphere represents a significant source of virus diversity and activity, indicating that viruses might have a considerable role in the ecology of the peat moss (Sphagnum) microbiome. Hühns et al. (2003) mechanically inoculated P. patens with tomato spotted wilt virus and detected viral proteins by Western blotting in P. patens gametophores. However, they could not detect symptoms such as chlorosis or wilting in P. patens gametophores. Polischuk, Budzanivska, Shevchenko, and Oliynik (2007) reported detecting antigens of tobacco mosaic virus and cucumber green mottle mosaic virus in samples of Barbilophozia and Polytrichum mosses, respectively.

Although these studies reveal genomic traces of past viral infections in bryophytes and highlight the natural omnipresence of viruses, the symptoms, interactions and developmental consequences of viral virulence factors, such as VSRs, remain unknown within this lineage. Here we characterized the phenotypic effects caused by two VSRs of vascular plants, 2b and P19, in P. patens. The analyses revealed that altered sRNA levels and accumulation of transcription factor mRNAs by VSR P19 has a substantial impact on the developmental processes of P. patens.

MATERIALS AND METHODS

Plant material and growth conditions

Protonemal tissue of Physcomitrella patens (Hedw.) Bruch & Schimp. (family Funariaceae, recently renamed to Physcomitrium patens (Hedw.) Mitt.) ecotype Gransden (Ashton & Cove, 1977; IMSC accession number 40001) was grown on cellophane‐covered BCD medium (1 mM MgSO4, 1.85 mM KH2PO4 [pH adjusted to 6.5 with KOH], 10 mM KNO3, 45 μM FeSO4, 0.22 μM CuSO4, 0.19 μM ZnSO4, 10 μM H3BO4, 0.10 μM Na2MoO4, 2 μM MnCl2, 0.23 μM CoCl2, 0.17 μM KI) (Ashton & Cove, 1977) supplemented with 1 mM CaCl2 and 45 μM Na2EDTA and solidified with 0.8% agar (Lehtonen et al., 2009) and on MM medium (BCD medium supplemented with 5 mM ammonium tartrate). Moss cultures were subcultured weekly. Protonemal tissue of P. patens was collected and homogenized in sterile Milli‐Q water (Millipore, Billerica, MA, USA) with an Ultra‐turrax T25 homogenizer (Janke & Kunkel IKA‐Labortechnik) and plated on cellophane‐covered BCD or MM medium in a growth chamber (Model 3,755; Forma Scientific, Marietta, OH, USA) at 20°C, with a 12 hr: 12 hr, light: dark photoperiod and light intensity of 40 μmol m−2 s−1.

Gateway cloning of inducible 2b and P19 suppressor constructs

Gateway cloning was used to obtain VSR P19 of tomato bushy stunt virus (TBSV) (gene ID, 1493957) and 2b of CMV (gene ID, 7552468; kin strain) (Siddiqui, Sarmiento, Turve, Lehto, & Lehto, 2008; Siddiqui, Valkonen, Rajamäki, & Lehto, 2011) under an estrogen inducible promoter in the pPGX8 targeting construct (Ishikawa et al., 2011). The CMV kin strain produces only a mild mosaic disease on tobacco whereas the Y strain causes severe yellow mosaics. For gateway entry clone preparation, attB1 and attB2 sequences were added to the 5′ end and to the 3′ end of the P19 and 2b genes by PCR amplification with attB1F and attB2R primers (10 μM each primer) (Table S1). PCR was carried out with the Phusion polymerase (Thermo Fisher Scientific, Vilnius, Lithuania) in a final volume of 50 μl. The PCR amplification products (303 bp for 2b and 520 bp for P19) were purified with an E.Z.N.A. Gel Extraction kit (Omega Bio‐tek, Inc., Georgia, USA) and quantified with a Nanodrop 2000 spectrophotometer (Thermo Scientific, Wilmington, USA).

Recombination of the PCR product and the Gateway donor vector pDONR221 (Thermo Fisher Scientific) was carried out in a reaction mixture containing 2 μl plasmid pDONR221 (75 ng μl−1), 2 μl of attB PCR product and BPII clonase (1 μl). Recombination was done at 25°C overnight, followed by Proteinase K treatment at 37°C for 10 min. Escherichia coli DH5α cells were transformed with the entry clones containing the P19 or 2b suppressor encoding region flanked by attB sequences. Plasmids were purified with the GenElute Plasmid Miniprep kit (Sigma, St. Louis, USA).

Entry clones (containing the gene of interest flanked by attL sites) were transferred to the gateway destination vector pPGX8 (containing attR sites) in the LR reaction. LR‐reactions refer to the recombination of entry clone (containing attL sites) and destination vector (containing attR sites) to generate an expression clone. The LR reaction mixture for each suppressor contained an entry clone (95 ng for 2b and 86 ng for P19), the destination vector pPGX8 (156 ng) and 1 μl of LRII clonase. The reaction mixture was incubated at 25°C for 1 hr, treated with Proteinase K for 10 min at 37°C and used to transform E. coli DH5α cells. pPGX8 containing the gene for viral suppressor 2b or P19 (the targeting constructs) was isolated from transformed E. coli DH5α with the Qiagen plasmid maxiprep kit (Qiagen, Hilden, Germany).

Transformation of

The targeting constructs were digested with PmeI (5 U/μl), and digestion products were purified by phenol/chloroform extraction. For transformation of P. patens, protoplasts were prepared according to Schaefer, Zryd, Knight, and Cove (1991). Protoplasts (300 μl of 1.2 × 106 protoplasts ml−1) were transformed with 10–30 μg of DNA and plated on cellophane‐covered BCD plates supplemented with 10 mM CaCl2, 45 μM Na2‐EDTA, 5 mM ammonium tartrate, 6.6% mannitol and 0.5% glucose. The growth conditions were the same as described for moss material above. After 6 days of incubation, protoplasts on the cellophane were transferred to MM‐medium containing hygromycin (30 μg ml−1; Fluka, Seelze, Germany). Selection was performed by two rounds of selection: regenerating protoplasts were grown for 2 week with hygromycin, followed by 2 week of growth without selection and then 2 week of growth with hygromycin.

The correct integration of the knock‐in construct was verified in the moss lines surviving the selection. Integration of the construct was verified by PCR with the PIG1 forward primer and the LexA reverse primer for the 5′ integration site and with the PIG1 reverse primer and 35S poly‐A forward primer (Table S1) for the 3′ integration. The Phire Plant Direct PCR Kit (Thermo Fisher Scientific) was used for screening mutant lines according to the manufacturer's instructions. Three P. patens P19 lines (P19‐14, P19‐26 and P19‐27) and three P. patens 2b lines (2b‐1, 2b‐10 and 2b‐21), whose phenotypes resembled wild type under non‐inducing conditions, were tested for ploidy levels according to Schween, Gorr, Hohe, and Reski (2003), and because of wild‐typical ploidy level, they were chosen for further experiments. In all experiments, wild‐type P. patens was included for comparison.

Chemical treatments

The estrogen receptor‐based chemical‐inducible system of Ishikawa et al. (2011) was used to test silencing suppression activity of the two VSRs, 2b and P19, in P. patens. A stock solution of β‐estradiol (10 mM; Sigma) was prepared in dimethyl sulfoxide (DMSO) and stored at −20°C until use. After 7 days of growth on cellophane‐covered MM plates, P. patens was treated with 10 ml of 1 μM β‐estradiol solution (10 mM stock diluted with sterile water) or with 10 ml of sterile water. Samples for RNA extraction were taken 7 and 12 days after treatment.

Protein extraction and tandem mass spectrometry (MS/MS)

Total protein was extracted according to Cove et al. (2009). In brief, P19‐14 and 2b‐1 moss tissue was ground in liquid nitrogen. The resulting powder was divided into 100‐μl aliquots, and 1 ml pre‐chilled protein extraction solution (0.07% [v/v] β‐mercaptoethanol, 1% [w/v] trichloroacetic acid in 96% ethanol) was added to each. The mixture was incubated at −20°C for 2 hr and centrifuged at 20,000g for 30 min at 4°C, followed by treatment with protein wash solution (2 mM EDTA, 0.07% [v/v] β‐mercaptoethanol, 100 mM phenylmethylsulfonyl fluoride, in acetone). Proteins were washed until the supernatant was completely clear and were then dissolved in 4 M urea.

Cysteine bonds were reduced with 0.05 M dithiothreitol for 20 min at 37°C and alkylated with 0.15 M iodoacetamide at room temperature. Samples were digested with 1 μg trypsin (sequencing grade modified trypsin, V5111, Promega, Madison WI, USA) overnight at 37°C, and peptides were then quenched with 10% trifluoroacetic acid (TFA) and purified with C18 microspin columns (Harvard Apparatus, USA) using 0.1% TFA in 50% acetonitrile to elute the peptides. The dried peptides were reconstituted in 30 μl of buffer A (0.1% TFA in 1% acetonitrile).

Liquid chromatography coupled with tandem mass spectrometry (LC–MS/MS) analysis was carried out on an EASY‐nLC1000 (Thermo Fisher Scientific) connected to a Velos Pro‐Orbitrap Elite hybrid mass spectrophotometer (Thermo Fisher Scientific) with a nano‐electrospray ion source (Thermo Fisher). Samples were separated using a two‐column setup.

The elution program was 5% buffer B for 5 min, then a linear gradient to 35% B over maintained for 60 min; then a linear gradient to 80% B over maintained for 5 min; and finally a linear gradient to 100% maintained for 10 min. The flow rate was 0.3 μl/min. A sample volume of 6 μl was injected per LC–MS/MS run. The MS2 method was used to fragment the 20 most intense precursor ions with collision‐induced dissociation (CID) (energy 35). Acquired MS2 scans were searched against a protein database (Varjosalo Laboratory) using the Sequest search algorithms in Thermo Proteome Discoverer (Thermo Fisher Scientific). The allowed mass error for the precursor ions was 15 ppm and was 0.8 Da for the fragments. A static residue modification parameter was set for the carbamidomethyl group +57.021 Da (C) for cysteine residues. Methionine oxidation was set as a dynamic modification of +15.995 Da (M). Only full tryptic peptides were considered (i.e., no fragments of tryptic peptides), and a maximum of one missed cleavage was allowed (false‐discovery rate 0.05).

RNA extraction

Total RNA was extracted according to Chang, Puryear, and Cairney (1993). Moss tissues were ground in liquid nitrogen and transferred to pre‐warmed isolation buffer (2% cetyl trimethyl ammonium bromide, 2% polyvinylpyrrolidone K‐30, 100 mM Tris–HCl [pH 8.0], 25 mM EDTA, 2 M NaCl, 0.02% β‐mercaptoethanol). Samples were extracted three times with chloroform/isoamyl alcohol (24:1), shaken at 250 rpm for 15 min and centrifuged at 10,000g for 15 min. Samples were precipitated with 8 M LiCl at 4°C overnight. After centrifugation at 10,000g for 30 min at 4°C, the pellet was dissolved in pre‐warmed SSTE buffer (1 M NaCl, 0.5% SDS, 10 mM Tris–HCl [pH 8.0], 1 mM EDTA) and extracted once with chloroform/isoamyl alcohol (24:1) as above. The RNA was precipitated from the aqueous phase with 2 vol absolute ethanol and dissolved in nuclease‐free water. The RNA concentration was measured with a spectrophotometer (Nanodrop 2000, Thermo Scientific).

Northern blotting

RNA (5 μg) extracted from wild‐type P. patens and the lines 2b‐1 and P19‐14 was analyzed by agarose gel electrophoresis on a 1% (w/v) formaldehyde gel, blotted onto Amersham Hybond nx nylon membrane (GE Healthcare) and crosslinked with UV light (Sambrook, Fritsch, & Maniatis, 2001). Plasmids containing the 2b or P19 suppressor sequences were used as DNA templates for probe amplification, which was done with Dynazyme II DNA polymerase (Thermo Fisher Scientific) with their respective primers (Table S1) according to the manufacturer's instructions. Probe sizes for 2b and P19 were 237 and 251 bp, respectively. The radioactive probes were prepared with 5 μl of template DNA, 13.5 μl nuclease‐free water and 0.5 μl random hexamers (200 ng/μl) and were boiled for 10 min, followed by cooling on ice for 5 min. Klenow buffer (3 μl), 1 μl Klenow fragment, 3 μl of dNTP‐dCTP and 4 μl of [α‐32P]dCTP (PerkinElmer, Turku, Finland) was added, followed by incubation at 37°C for 1 hr. The probe was purified with the QIAquick Nucleotide Removal kit (Qiagen). Hybridization was carried out at 65°C overnight. The membrane was washed three times with 1× SSC (3.0 M NaCl and 0.3 M sodium citrate), placed into a cassette with an imaging plate and exposed overnight. Radiation energy was scanned and documented with the fluorescent image analyzer FLA‐5100 (Fuji Photo Film Co. Ltd., Tokyo, Japan).

mRNA sequencing

Twelve RNA sequencing libraries were created in total: samples were taken at 7 and 12 days after β‐estradiol or water treatments from the suppressor lines 2b‐1 and P19‐14 as well as from wild‐type P. patens. For RNA sequencing, 1 μg of RNA was treated with Ribo‐Zero (Illumina) to remove ribosomal RNA. The NEBNext Ultra Directional RNA Library Prep Kit (Illumina) was used to generate a cDNA library. Briefly, the adapters for sequencing (Illumina) were ligated, and cDNA of the correct size was selected with AMPure XP Beads. USER Enzyme (uracil‐specific excision reagent) was used to remove deoxyuridine residues. The cDNA library was amplified with high‐fidelity PCR and purified with AMPure XP Beads. The quality of the library was assayed with an Agilent 2100 Bioanalyzer (Agilent) and library quantity was determined by Qubit (Qubit Fluorometer V.1.27). Paired‐end (2 × 75 bp) RNA sequencing was performed at the Functional Genomics Unit (Helsinki) via the Illumina NextSeq sequencing system.

Processing of RNA sequences and identification of candidates for qPCR

Adapter removal was performed as part of the sequencing process by the service provider, Functional Genomics Unit. Paired‐end reads were mapped to the P. patens genome V3 (Lang et al., 2018) with STAR version 2.6.0c (Dobin et al., 2013) and only concordantly mapped read pairs were retained. htSeqCount version 0.11.0 (Anders, Pyl, & Huber, 2015) was then used in union mode to count the reads that mapped to exons of P. patens V3.3 gene models (Lang et al., 2018) while ignoring secondary or supplementary alignments. The derived feature counts served as input for an explorative gene expression analysis based on R and the DESeq2 package (R Core Team, 2021; Love, Huber, & Anders, 2014) to identify candidate genes for quantitative RT‐PCR.

cDNA synthesis

RNA samples from two biological replicates were pooled, and total RNA was treated with DNase (1 U/μl, Promega). For cDNA synthesis, RNA (1 μg) treated with DNase was incubated with 1 μl (200 ng/μl) random hexamers (Thermo Fisher Scientific) at 65°C for 5 min. Reverse transcription was performed using 4 μl of 5× reaction buffer (250 mM Tris–HCl [pH 8.3 at 25°C], 250 mM KCl, 20 mM MgCl2, 50 mM DTT [Thermo Fisher Scientific]), 0.5 μl riboblock RNase inhibitor (40 U/μl, Thermo Fisher Scientific), 2 μl dNTP[s] (10 mM each, Thermo Fisher Scientific) and 1 μl RevertAid M‐MuLV reverse transcriptase (200 U, Thermo Fisher Scientific). Reaction mixtures were incubated at 25°C for 10 min, followed by incubation at 42°C for 1 hr. The reaction was stopped by heating at 70°C for 10 min, and the resulting cDNA was diluted with 80 μl of nuclease‐free water.

Quantitative RT‐PCR

For primer design, selected P. patens V3.3 gene sequences were obtained from Phytozome (Goodstein et al., 2012; Lang et al., 2018). Primers for qPCR were designed with GenScript (www.genscript.com/tools/pcr-primers-designer). Target specificities of the primers were confirmed by sequencing the amplification product. The gene encoding ribosomal subunit L21 was used as a reference gene (Lehtonen, Akita, Frank, Reski, & Valkonen, 2012) to normalize variation in cDNA amounts. Quantitative PCR reactions were carried out using a reaction volume of 10 μl, which included 0.5 μl of the forward and reverse primers (5 μM each), 4 μl of cDNA sample and 5 μl of LightCycler 480 SYBR Green I Master mix (Roche Diagnostics GmbH Germany). Three technical replicates of each sample were included. Quantitative PCR was carried out according to the manufacturer's instructions. The relative expression ratio of each target gene was calculated according to Pfaffl (2001). Statistical analyses were carried out to define whether there were significant differences in the expression levels among 2b and P19 suppressor lines and wild‐type P. patens. Fold values were log2‐transformed and then subjected to ANOVA analysis (IBM SPSS, v. 25.0, SPSS Inc.).

sRNA extraction and sequencing

For sRNA extraction, wild‐type P. patens and lines P19‐14 and 2b‐1 were grown for 7 days on cellophane‐covered MM plates and treated with β‐estradiol. Samples were taken 7 days after β‐estradiol treatment, frozen in liquid nitrogen and stored at −80°C prior to extraction. sRNAs were extracted by using the miRNeasy Mini kit (miRNeasy Mini Handbook 12/2014, Qiagen, GmbH Hilden, Germany), and DNase treatment was performed with the RNase‐Free DNase set (Qiagen). Extracted sRNAs were eluted into 30 μl of nuclease‐free water. The experiment was replicated three times. RNA samples were sent to Fasteris SA (Plan‐les‐Ouates, Switzerland) for sequencing on an Illumina HiSeq 3000/4000 system.

sRNA data processing and miRNA expression analysis

The sRNA sequencing libraries comprised 25,508,567 (WT replicate 1), 31,285,278 (WT replicate 2), 47,707,278 (WT replicate 3), 33,434,271 (2b replicate 1), 41,837,206 (2b replicate 2), 46,889,931 (2b replicate 3), 42,172,120 (P19 replicate 1), 37,003,119 (P19 replicate 2) and 41,335,801 (P19 replicate 3) single‐end raw reads. Prior to mapping, they were adapter‐trimmed and freed of empty inserts using Trimmomatic (Version 0.32; Bolger, Lohse, & Usadel, 2014). sRNA clusters were identified with ShortStack (Version 3.8.5; Johnson, Yeoh, Coruh, & Axtell, 2016), which allowed for a single mismatch of reads mapped to the P. patens genome with bowtie (Langmead, Trapnell, Pop, & Salzberg, 2009) and used the default ‘unique – seeded guide’ method for automatic placement of multi‐mapping reads. Reference miRNA precursors from P. patens for cluster annotation were obtained from miRbase (Release 22.1; Kozomara, Birgaoanu, & Griffiths‐Jones, 2019). The output of ShortStack was parsed to retrieve a consensus set of unique, non‐overlapping clusters that passed ShortStack's internal miRNA classification. Based on the mapped read counts reported by ShortStack we performed a differential expression analysis of this consensus set using R and the DESeq2 package (Love et al., 2014; R Core Team, 2021). Expressed sequences of the three contrasts P19 vs. WT, P19 vs. 2b and 2b vs. WT were filtered for an adjusted p‐value of <0.05 and a log2 fold‐change of <−0.58 or >0.58, respectively.

miRNA target prediction and gene ontology enrichment

Identified differentially expressed miRNAs were passed on to three target prediction tools, psRNAtarget (Dai & Zhao, 2011), TAPIRhybrid (Bonnet, He, Billiau, & van de Peer, 2010) and TargetFinder (Fahlgren et al., 2007) with default parameters. The output was parsed and subsequently filtered to obtain miRNA‐target‐pairs found by at least two predictors. These target transcripts and their corresponding annotation (Lang et al., 2018) were then subject of a gene ontology (GO; Ashburner et al., 2000; The Gene Ontology Consortium, 2019) enrichment analysis based on R and the topGO package (Alexa, Rahnenführer, & Lengauer, 2006; R Core Team, 2021). The top 25 enriched terms for each comparison and each GO category were obtained by minimum weight as the result of the built‐in Fisher test. Overlaps of miRNA or targets and word clouds of enriched terms were plotted using the packages upSetR (Conway, Lex, & Gehlenborg, 2017) and GO summaries (Kolde & Vilo, 2015).

RESULTS

Induced expression of P19 and 2b suppress development of leafy gametophores

Transgenic lines of P. patens with estrogen receptor‐based β‐estradiol‐inducible VSR P19 of TBSV or 2b of CMV were generated using the pPGX8 plasmid. Integration of the constructs was verified with PCR using the PIG1 (forward) and LexA (reverse) primers for the 5′ integration site and the PIG1 (reverse) and 35SpolyA (forward) primers for the 3′ integration site (Table S1, Figure S1). Three P19 suppressor lines (P19‐14, P19‐26 and P19‐27) and three 2b suppressor lines (2b‐1, 2b‐10 and 2b‐21) had ploidy levels that were the same as wild type and thus they were used for further analyses.

Water‐treated and β‐estradiol‐treated wild‐type P. patens developed normally, and adult leafy gametophores were observed 12 days after water treatment (Figure 1). Water‐treated P19 (P19‐14, P19‐26 and P19‐27) and 2b (2b‐1, 2b‐10 and 2b‐21) suppressor lines also developed normally, similar to wild‐type P. patens (Figure 1). However, the formation of leafy gametophores was completely disrupted in suppressor lines expressing the P19 protein and was reduced to a lesser extent in the 2b suppressor line, but not in wild‐type P. patens indicating that P19 and 2b interfered with the development of P. patens.

Figure 1

Phenotype of Physcomitrella patens 2b and P19 suppressor lines after water and β‐estradiol treatments. (a) Overview of plants treated with water or β‐estradiol 14 days after treatments. Wild‐type (WT) P. patens and 2b‐1 and P19‐14 suppressor lines developed leafy gametophores 14 days after water treatment. The development of leafy gametophores in β‐estradiol‐treated suppressor lines was clearly reduced: the P19‐14 suppressor line did not develop leafy gametophores, and in the 2b‐1 suppressor line the formation of leafy gametophores was reduced and they were smaller when compared with wild‐type P. patens. (b) Microscopic observation revealed that 7 days after water treatment all plants had developed buds (black arrows). Buds were also observed in β‐estradiol‐treated wild type and the 2b‐1 suppressor line, but distinct buds or leaf initials were not observed in the P19‐14 suppressor line. Furthermore, the β‐estradiol‐treated P19‐14 line had shorter protonema cells than water‐treated plants. All plants developed leafy gametophores by 12 days after water treatment. Wild‐type plants had also developed normally by 12 days after β‐estradiol treatment. The development of gametophores was delayed in the β‐estradiol‐treated 2b‐1 line and was completely impaired in the P19‐14 suppressor line. Scale bar 100 μm [Colour figure can be viewed at wileyonlinelibrary.com]

A Northern blotting experiment confirmed the absence of transcripts for P19 and 2b suppressors in water‐treated wild‐type P. patens, the P19‐14 and 2b‐1 lines and also in β‐estradiol‐treated wild‐type plants (Figure 2). We did, however, observe transcripts for P19 and 2b in β‐estradiol‐treated P19‐14 and 2b‐1 lines 7 days after treatment (Figure 2), confirming the expression of P19 and 2b following β‐estradiol treatment. Furthermore, the protein sequence of the P19 suppressor was detected by tandem mass spectrometry after β‐estradiol‐treatment, indicating the presence of P19 suppressor at the protein level in P. patens line P19‐14 (Table S2). We could not detect the protein sequence of the 2b suppressor in P. patens line 2b‐1 by tandem mass spectrometry.

Figure 2

Northern blot analysis of Physcomitrella patens lines transformed with vectors containing viral suppressor 2b or P19. Foreign transcripts are inducible by β‐estradiol in this system. Plants treated with β‐estradiol accumulated 2b and P19 transcripts (black arrows), whereas plants treated with water did not show such expression. Wild‐type plants treated with β‐estradiol did not show any accumulation of 2b‐1 or P19‐14 transcripts or unspecific hybridization [Colour figure can be viewed at wileyonlinelibrary.com]

Gene expression analysis reveals candidate genes affected by viral suppressors

Unreplicated RNA‐seq analysis was performed for screening the candidate genes for qPCR experiment. Single RNA expression libraries representing treated and untreated transformed lines or wild‐type at two separate time points (12 in total) were created and sequenced for discovering VSR target genes for subsequent analyses. RNA sequencing reads were mapped against the P. patens genome with efficiencies between 85.7% and 89.7% resulting in 21.4 million to 31.8 million uniquely mapped reads for the individual libraries. Based on these mappings we identified 20,043–21,115 genes that had read counts in the individual samples, and variation in their expression was evaluated via a principal component analysis (Figure S2), revealing distinct differences in β‐estradiol‐induced P19‐ and 2b‐expressing lines compared with controls of untreated transformed or wild‐type lines, respectively. P19 samples that had been treated for 7 and 12 days shared high similarity and were clearly separated from the two clusters comprised of the individual controls for the two distinct time points. Concerning 2b, variation was highest between the individual time points, whereas only the β‐estradiol‐induced sample taken after 12 days was distinctly separate from all other samples. Gene expression levels were compared between each of the induced samples and its three respective controls. We then calculated the mean log2‐transformed fold change for genes that showed a consistent direction of regulation in all three individual comparisons for each induced sample. Subsequently, we compiled consensus lists of genes exhibiting highly altered expression levels by retaining only those genes with a mean log2 fold change of either >5 or <−5, likely representing targets of viral suppressors (Table S3).

Silencing suppressors interfere with the expression of transcription factor genes

We aimed to find candidate genes that have been described to play a role during moss development, are putative targets of miRNAs and whose expression level was strongly altered upon expression of P19 or 2b (Table S4) for further validation via qPCR. By combining the results from our RNA‐seq analysis with functional gene annotations (Lang et al., 2018) and gene expression data from the Physcomitrella gene atlas (Perroud et al., 2018), we identified four transcription factors as candidates and were able to confirm their altered expression levels: Pp3c13_3710V3.1 is an AP2/erf domain‐containing transcription factor, Pp3c1_40330V3.1 is the ethylene‐responsive transcription factor EREB‐11 and Pp3c1_42920V3.1 and Pp3c23_5390V3.1 are predicted MYB transcription factors. Expression levels of the four transcription factors were tested by qRT‐PCR in P19 (P19‐14, P19‐26 and P19‐27) and 2b (2b‐1, 2b‐10 and 2b‐21) suppressor lines and also in wild‐type P. patens (Figure 3).

Figure 3

Transcript levels of four transcription factors were determined in suppressor lines P19 (P19‐14, P19‐26 and P19‐27) and 2b (2b‐1, 2b‐10 and 2b‐21) and in wild‐type Physcomitrella patens. The transcript levels of four transcription factors, (a) Pp3c13_3710V3.1, (b) Pp3c1_40330V3.1, (c) Pp3c1_42920V3.1 and (d) Pp3c23_5390V3.1, were quantified across three experiments. Statistical analyses (ANOVA, Tukey's test) showed that expression of all four transcription factor genes was significantly higher in β‐estradiol‐treated P19 suppressor lines than in β‐estradiol‐treated wild‐type plants or β‐estradiol‐treated 2b suppressor lines (p‐value <0.05). Black bars indicate β‐estradiol treatment, and gray bars indicate water treatment. Error bars indicate the log2‐transformed SD of fold change from maximum value. X‐axis labels indicate the moss line, treatment and time in which samples were taken (mq: milliQ water; 7 days: samples were taken 7 days after treatment; 12 days: samples were taken 12 days after treatment)

We found no statistically significant differences in the expression levels of the four candidate transcription factor genes when comparing water‐treated wild‐type plants with P19 or 2b suppressor lines. In water‐treated samples, the only statistically significant difference was observed between water‐treated P19 and 2b lines for the Pp3c1_42920V3.1 gene. The transcript levels of the four transcription factors were clearly higher in P19 suppressor lines treated with β‐estradiol (Figure 3) compared with water‐treated P19 suppressor lines, water‐treated or β‐estradiol‐treated wild‐type plants or 2b suppressor lines. Statistical analyses (ANOVA, Tukey's test) showed that gene expression levels of all four tested transcription factors were indeed significantly higher in β‐estradiol‐treated P19 suppressor lines than in β‐estradiol‐treated wild‐type plants or β‐estradiol‐treated 2b suppressor lines (p < 0.05) (Table 1).

Table 1

Statistical analyses of the fold change in transcription factor transcripts

Transcription factor	Transcription factor family	ANOVA
		Water treatment		β‐estradiol treatment
		F‐value	Sig.	F‐value	Sig.
Pp3c1_42920V3.1	MYB	3.69	0.033	7.22	0.002
Pp3c23_5390V3.1	MYB	3.02	0.059	6.30	0.004
Pp3c1_40330V3.1	AP2/EREBP	0.57	0.569	25.2	0.000
Pp3c13_3710V3.1	AP2/EREBP	2.30	0.112	7.20	0.002

P19 expression alters miRNA levels in moss

As P19 and 2b affect RNA interference in vascular plants during viral infection, we tested for differences in the sRNA pool of P. patens upon their synthesis in this moss. Therefore, we sequenced sRNAs from wild type as well as from P19‐ and 2b‐expressing lines, each in triplicate. Examination of the length profiles of reads mapped to the P. patens genome revealed a strong bias toward 20/21 nucleotides for the P19 line compared with 2b and wild type (Figure 4a). We therefore binned the mapped reads into four categories based on their lengths: CatI (<20 nt), CatII (20–22 nt), CatIII (23–25 nt) and CatIV (>25 nt) (Figure 5). A two‐factor ANOVA then confirmed statistically significant differences between P19 and wild type (p < 0.05), most prominently between P19‐CatII and WT‐CatII (p < 1E‐07), but not between 2b and wild type.

Figure 4

sRNA sequencing reveals a strong bias in miRNAs in the P19‐expressing line of Physcomitrella patens compared to wild type (WT). (a) Principal component (PC) analysis of mapped read lengths shows a bias for 20/21 nucleotides for the P19 replicates. (b) PC analysis of differential miRNA expression. (c) UpSet plot to highlight the overlap of up‐ and down‐regulated miRNAs in the three comparisons P19 vs. WT, P19 vs. 2b and 2b vs. WT. (d) UpSet plot of predicted targets of differentially expressed miRNAs of (c) [Colour figure can be viewed at wileyonlinelibrary.com]

Figure 5

Distribution of mapped sRNA read lengths. The P19‐expressing line reveals a strong shift towards reads between 20 and 22 nt in length, whereas reads of lengths shorter or longer are underrepresented in comparison to the 2b‐expressing line and wild type (WT). n = 3 [Colour figure can be viewed at wileyonlinelibrary.com]

To investigate the effects on sRNA composition and on the abundance of miRNAs in particular, we performed two independent runs of sRNA gene annotations with ShortStack (Johnson et al., 2016): a reference‐free de novo identification of sRNA clusters and a reference‐based miRNA annotation by supplying ShortStack with all 262 primary miRNA sequences of P. patens publicly available in miRBase (Kozomara et al., 2019). By merging both predictions we obtained a set of 139 unique miRNA precursor genes (Supplemental File 1 in Appendix S1). In alignment with our findings on sRNA length composition in the three lines, a subsequent differential expression analysis revealed that the P19 line clearly separated from the other two lines in its individual expression pattern of miRNAs (Figure 4b). We found 92 and 83 miRNAs to be significantly differentially expressed in the comparisons of P19 vs. WT and P19 vs. 2b, respectively, while we found only 6 when comparing 2b with WT (Supplemental File 1 in Appendix S1). Among the former we found homologs of ppt‐miR160 and ppt‐miR166, two widely conserved plant miRNA families, as well as members of the four moss‐specific families ppt‐miR534, ppt‐miR1023, ppt‐miR898 and ppt‐miR902. While ppt‐miR160b, ppt‐miR160d, ppt‐miR160h, ppt‐miR160i, ppt‐miR166c, ppt‐miR166d, ppt‐miR166f, ppt‐miR166i, ppt‐miR166l, ppt‐miR534a, ppt‐miR534b, ppt‐miR898a, ppt‐miR1023a, ppt‐miR1023b, ppt‐miR902d, ppt‐miR902f, ppt‐miR902g were downregulated in the P19‐expressing line, ppt‐miR1023c and ppt‐miR902j were upregulated. The putative targets of these miRNAs have been described previously as an auxin response factor in the case of miR160, a class III HD‐ZIP gene for miR166, a BLADE ON PETIOLE 2‐like BTB gene and ankyrin‐domain genes for miR534, K‐box containing transcription factors for miR898, a protein kinase for miR902 and the target of miR1023 has been described to be similar to STOP1 (sensitive to proton rhizotoxicity) (Addo‐Quaye et al., 2009; Arazi, 2012; Axtell, Snyder, & Bartel, 2007; Floyd & Bowman, 2004).

Targets of altered miRNAs in P19 are involved in fundamental processes

To investigate the response to differential miRNA expression, we obtained a high‐confidence set of putative miRNA targets by performing three independent predictions with psRNAtarget (Dai & Zhao, 2011), TAPIRhybrid (Bonnet et al., 2010) and TargetFinder (Fahlgren et al., 2007) as their prediction overlap was reported to be reliable (Srivastava, Moturu, Pandey, Baldwin, & Pandey, 2014). Corresponding to 84, 75 and all 6 differentially expressed miRNAs in the comparisons between P19 and WT, P19 and 2b as well as 2b and WT, we were able to predict 362, 308 and 39 target transcripts encoded by the P. patens genome (Supplemental File 1 in Appendix S1). Just as for the differentially expressed miRNAs themselves (Figure S1c), we found the biggest overlap of predicted targets to be between comparisons involving P19, in both directions of gene expression (Figure S1d), whereas comparisons between 2b and WT revealed only few targets.

In order to gain insights into the consequences of altered miRNA expression in our lines, we performed a gene ontology (GO; Ashburner et al., 2000, The Gene Ontology Consortium, 2019) enrichment analysis of the predicted targets. Therefore, we obtained GO annotations for the target transcripts from Lang et al. (2018) and identified the top 25 of enriched nodes for each of the three comparisons and the three GO categories ‘biological process’ (BP), ‘molecular function’ (MF) and ‘cellular component’ (CC; Supplemental File 1 in Appendix S1). In the comparisons with both 2b and WT, the top enriched terms corresponding to targets of up‐regulated miRNAs in the P19 line comprise among others protein phosphorylation, metabolic and biosynthetic processes or leaf development, whereas terms relating to transcriptional regulation, auxin signalling or light detection and phototransduction are among those enriched in targets of down‐regulated miRNAs (Figure 6, Figures S3 and S4).

Figure 6

Gene ontology (GO) enrichment identifies miRNA targets are involved in metabolic, signalling and developmental processes. Word clouds of enriched GO terms in the category ‘biological process’ of transcripts targeted by differentially expressed miRNAs. Coloring corresponds to the terms' weight of topGO Fisher Test and size corresponds to individual rank in each enrichment

DISCUSSION

Our present study reveals that the expression of VSR P19 of TBSV alters the development, gene expression and the sRNA pool of the moss P. patens. We found that four different transcription factors (AP2/erf, EREB‐11 and two MYBs) accumulated in the P19 lines. sRNA sequencing revealed that VSR P19 significantly changed the microRNA pool in P. patens. However, we could not find significant differences between wild type and transgenic 2b line in qPCR experiments. Northern blot and mass spectrometry analysis confirmed the expression of P19 at transcript and protein level. P. patens plants transformed with VSR 2b of CMV accumulated 2b transcript, but we could not detect 2b protein. Failing protein detection following heterologous expression in moss has recently been observed by Top et al. (2021), who described the phenomenon as ‘heterosplicing’. Heterosplicing, which may result in accumulation of transcripts but not in functional protein in P. patens, could be a reason why we were unable to detect 2b protein in our current experiments, although we applied highly sensitive mass spectrometry.

P. patens plants transformed with VSR P19 of TBSV were developmentally impaired when compared to wild‐type plants as they did not develop leafy gametophores. Transgenic lines of the vascular plant N. benthamiana expressing viral RNA silencing suppressor P19 of TBSV showed blistering or necrotic lesions within the leaves (H. B. Scholthof, Scholthof, & Jackson, 1995; Siddiqui et al., 2008), whereas N. benthamiana plants expressing the P19 gene of cymbidium ringspot virus have developmental disorders with leaf distortions and elongated stem internodes (Kontra et al., 2016). These results indicate the VSR P19 induces developmental disorders both in vascular and non‐vascular plants.

Studies on VSR‐expressing transgenic plants have shown that viral suppressors strongly interfere with endogenous silencing pathways (Chapman, Prokhnevsky, Gopinath, Dolja, & Carrington, 2004; Kasschau et al., 2003). Recently, Kontra et al. (2016) and Pertermann et al. (2018) reported that P19 efficiently binds endogenous miRNA duplexes, followed by the elevation of most of the mRNA targets as a consequence of miRNA duplex sequestration by P19 and the inability to program RISC for cleavage. Our qPCR and RNA‐seq data clearly showed that P19 of TBSV caused significant accumulation of transcripts of four different transcription factors belonging to MYB and APETALA 2/ethylene‐responsive element binding factor (AP2/EREBP) gene families in P. patens. In plants, both AP2/EREBP and MYB transcription factors regulate fundamental processes such as development, differentiation, metabolism and defense against pathogens (Ambawat, Sharma, Yadav, & Yadav, 2013; Dietz, Vogel, & Viehhauser, 2010). In addition, AP2/EREBP transcription factors are also important in the stress response in P. patens (Hiss et al., 2014; Richardt et al., 2010). Richardt et al. (2010) found the highest rates of induction under salt stress and/or abscisic acid for genes belonging to, among others, the AP2/EREBP family, whereas Jofuku, den Boer, Van Montagu, and Okamuroa (1994) found that the AP2 gene is expressed in the inflorescence meristem and in non‐floral organs including leaves and stem in Arabidopsis suggesting its crucial role in Arabidopsis development. Leech, Kammere, Cove, Martin, and Wang (1993) showed that in P. patens the corresponding transcript levels of two MYB‐like proteins, Pp1 and Pp2, reached maximum levels in young protonemal tissue that correlated with the time of maximum mitotic index. They also found that 1‐naphthaleneacetic acid resistant (nar) mutant lines, which are blocked at specific stages of development and are unable to produce gametophores, had aberrant levels of Pp1 and Pp2 transcripts, indicating that the normal expression of Pp1‐ and Pp2‐encoding genes is essential for cell growth during gametophytic development in P. patens. Our study clearly shows that P19 of TBSV caused significant accumulation of transcripts of MYB and AP2/EREBP gene families, indicating that the altered expression levels of these transcription factors might be linked to the developmental disorder observed in P. patens.

In plants, miRNAs are important regulators of multiple developmental processes (Axtell et al., 2007; Axtell & Bartel, 2005; Fattash, Voss, Reski, Hess, & Frank, 2007; Floyd & Bowman, 2004; Khraiwesh et al., 2010; Kindner & Martienssen, 2004; P. P. Liu et al., 2007; Mallory, Bartel, & Bartel, 2005). We found that P19 VSR alters the sRNA repertoire in the moss Physcomitrella. sRNA sequencing revealed a strong bias towards reads with a length of 20/21 nucleotides for the P19 line. sRNA predictions then verified significant differences in miRNA expression within this line with a total of 92 and 83 differentially expressed sequences in comparison to 2b and wild type, respectively. Interestingly, among miRNAs down‐regulated in the P19 line, we found the conserved sequences miR160a and members of the miR166 family indicating an association of conserved miRNAs with P19 or its targets, which has been reported previously in vascular plants (Kontra et al., 2016).

Our sRNA sequencing data also revealed that ppt‐miR534a, ppt‐miR534b and ppt‐miR1023a are down‐regulated in P19‐expressing P. patens when compared to wild type. Both ppt‐miR534a and ppt‐miR1023a are involved in the regulation of developmental processes of P. patens (Arazi, 2012). In Physcomitrella, miR534 targets mRNAs of blade on petiole (BOP) like proteins (Addo‐Quaye et al., 2009; Saleh et al., 2011); BOPs are transcriptional activators involved in cell differentiation in Arabidopsis and Physcomitrella (Ha et al., 2003; Hepworth, Zhang, McKim, Li, & Haughn, 2005; Jun, Ha, & Fletcher, 2010; Saleh et al., 2011). ppt‐miR1023a is highly abundant in 14‐days‐old protonema with young gametophores, whereas its abundance is much lower in young protonema and mature gametophores, indicating its specificity for buds (Arazi, 2012). The P19‐expressing P. patens line had severe developmental malformations. As ppt‐miR534a and ppt‐miR1023a are involved in cell differentiation in P. patens, down‐regulation of ppt‐miR534a and ppt‐miR1023a might explain the developmental arrest of the P19‐expressing P. patens line.

Our results showed that the level of miRNAs such as ppt‐miR160, ppt‐miR166, ppt‐miR534 and ppt‐miR1023 and the expression level of four transcription factors belonging to AP2/EREBP and MYB families were significantly altered in the P19‐expressing P. patens line. This generally indicates that the expression of the viral P19‐encoding region causes changes at both the miRNA and mRNA levels. miR160, miR166, miR534 and miR1023, as well as AP2/EREBP and MYB, regulate developmental processes in plants. We assume that P19 in particular causes the observed developmental disorders in the moss P. patens by altering key small RNAs and thereby modulating expression of genes involved in developmental processes.

P. patens has become a popular model system over the last decades, because of the many advantages it offers: axenic cultures are easy and rather fast to grow in Petri dishes, highly effective homologous recombination enables generation of gene replacement lines and microscopic examination of single cell‐layered tissues in the leafy gametophyte provides excellent possibilities for visualization (Reski, 2018). In addition to evolutionary studies and advances in understanding plant physiology and molecular biology using P. patens, numerous studies have proven the moss to be a suitable model also for phytopathology, including studies with fungi (Lehtonen, Akita, et al., 2012; Lehtonen et al., 2009; Lehtonen, Marttinen, Akita, & Valkonen, 2012; Ponce de León et al., 2007), oomycetes (Oliver et al., 2009; Ponce De León et al., 2012; Resemann et al., 2021) and bacteria (Andersson, Akita, Pirhonen, Gammelgård, & Valkonen, 2005; Ponce de León et al., 2007). RNA silencing pathways of bryophytes and vascular plants have significant similarities (You et al., 2017) – the current results indicate that the advantages of the Physcomitrella model can be utilized also in the field of plant–virus interaction studies.

CONFLICT OF INTEREST

The authors declare no conflicts of interest.

AUTHOR CONTRIBUTIONS

Mikko T. Lehtonen: Conceptualization; investigation; validation; visualization; writing‐review & editing. Eeva M. Marttinen: conceptualization; investigation; validation; visualization; writing mainly original draft, review & editing. Nico van Gessel: Conceptualization; data curation; formal analyses; visualization; writing partly original draft, review & editing. Ralf Reski: Conceptualization; funding acquisition; supervision; writing‐review & editing. Jari P. T. Valkonen: Conceptualization; resources; funding acquisition; supervision; writing‐ review & editing.

Figure S1. Verification of 2b and P19 transgenic moss lines by PCR. (a) Agarose gel analysis of 5′ integration site amplified with PIG1 fwd and LexA rev primer pair. (b) Agarose gel analysis of 3′ integration site amplified with 35 PolyA fwd and PIG1 rev primer pair. Wild‐type plants were used as negative control.

Click here for additional data file.

Figure S2. Principal component (PC) analysis reveals greater variability in gene expression in the P19 line of Physcomitrella patens than in the 2b line upon induction of viral suppressors. (A) Clear separation of the two β‐estradiol‐treated P19 samples, harvested 7 days (P19_E_7) and 12 days (P19_E_12) after treatment, from clusters of the two control groups comprising untreated P19 samples as well as treated (E) and untreated (C) wild‐type (WT) samples, harvested 7 days (WT_C_7, WT_E_7, P19_C_7) and 12 days (WT_C_12, WT_E_12, P19_C_12) after treatment. (B) Only the β‐estradiol‐treated 2b sample harvested after 12 days (2b_E_12) shows clear expression variance from either of the controls, whereas the other treated 2b sample (2b_E_7) is in a cluster with the control samples harvested after 7 days (WT_C_7, WT_E_7, 2b_C_7). Variability in expression mainly led to two distinct clusters comprising samples from the two harvesting time points.

Click here for additional data file.

Figure S3. Gene Ontology (GO) enrichment identifies miRNA targets are involved in molecular function. Colouring corresponds to the terms' weight of topGO Fisher test and size corresponds to individual rank in each enrichment.

Click here for additional data file.

Figure S4. Gene Ontology (GO) enrichment identifies miRNA targets are involved in cellular components. Coloring corresponds to the terms' weight of topGO Fisher test and size corresponds to individual rank in each enrichment.

Click here for additional data file.

Appendix S1. Supporting information.

Click here for additional data file.

Table S1. Primers used in this study.

Table S2. Protein data.

Table S3. DEGS.

Table S4. Fold change of genes that are potential targets of miRNAs.

Table S5. A complete list of available sequencing data sets.

Click here for additional data file.