Alstroemeria yellow spot virus (AYSV): a new orthotospovirus species within a growing Eurasian clade

An orthotospovirus distinct from all other orthotospoviruses was isolated from naturally infected alstroemeria plants. Disease symptoms caused by this virus mainly consisted of yellow spots on the leaves based on which the name alstroemeria yellow spot virus (AYSV) was coined. A host range analysis was performed and a polyclonal antiserum was produced against purified AYSV ribonucleoproteins which only reacted with the homologous antigen and not with any other (established or tentative) orthotospovirus from a selection of American and Asian species. Upon thrips transmission assays the virus was successfully transmitted by a population of Thrips tabaci. The entire nucleotide sequence of the M and S RNA segments was elucidated by a conventional cloning and sequencing strategy, and contained 4797 respectively 2734 nucleotides (nt). Simultaneously, a next generation sequencing (NGS) approach (RNAseq) was employed and generated contigs covering the entire viral tripartite RNA genome. In addition to the M and S RNA nucleotide sequences, the L RNA (8865 nt) was obtained. The nucleocapsid (N) gene encoded by the S RNA of this virus consisted of 819 nucleotides with a deduced N protein of 272 amino acids and by comparative sequence alignments to other established orthotospovirus species showed highest homology (69.5% identity) to the N protein of polygonum ringspot virus. The data altogether support the proposal of AYSV as a new orthotospovirus species within a growing clade of orthotospoviruses that seem to share the Middle East basin as a region of origin.

Introduction
Members of the Tospoviridae continuously cause significant losses in both mono-and dicot crops worldwide [1,2]. Although exact numbers indicating the economic impact of orthotospoviruses are not available, annual yield losses have been estimated to exceed at least 1 billion US$ for the type species Tomato spotted wilt virus (TSWV) only [3]. Due to this TSWV ranks high among plant viruses of economic importance worldwide [1,[4][5][6]. Orthotospoviruses are widespread and transmitted by at least 14 different thrips species [2,7] in a circulative and propagative manner. Virus particles are membrane-bound with a quasi-spherical morphology of about 80-120 nm in diameter and contain a tripartite ssRNA genome of negative and ambisense polarity [6]. The L RNA is of entire negative polarity and codes for only one gene in the viral complementary (vc) sense, the RNA-dependent RNA polymerase (RdRp; L protein) required for viral replication. The M and S RNA both have an ambisense gene arrangement and contain two ORFs separated by a non-coding intergenic region (IGR). The M RNA encodes a non-structural (NS M ) protein, in the viral (v) sense, that is involved in cell-to-cell movement and a glycoprotein precursor (G N -G C ) in the vc-sense required for (receptor-mediated) entry into the midgut-epithelium of thrips and vector-specificity [8][9][10][11][12][13]. The S RNA similarly encodes the non-structural protein (NSs) in the v-sense that Handling Editor: F. Murilo Zerbini.

Electronic supplementary material
The online version of this article (https ://doi.org/10.1007/s0070 5-018-4027-z) contains supplementary material, which is available to authorized users. 1 3 acts as suppressor of RNA silencing and supports efficient virus transmission by the vector, and the major structural nucleocapsid protein (N) in vc-sense [14][15][16][17][18]. Although host range, vector specificity and serology of the N protein are often used to identify and/or characterize orthotospoviruses, they only become assigned as a new species if their N protein exhibits less than 90% sequence homology to any of the established orthotospovirus species [2,19,20].
Ever since TSWV re-emerged in the 1980's and started its world wide spread by the introduction and spread of Frankliniella occidentalis, new viral species have been discovered. This has led to the identification of a current total of 29 orthotospovirus species of which 11 are officially recognized by the ICTV while the remaining, listed as tentative species [2,[20][21][22][23][24][25][26]. During the past two decades many viruses have been observed and identified in ornamental crops of which some cause significant economical losses [27].
In 2000 an alstroemeria plant showing only chlorotic and yellow spots on the leaves (Figure 1), and reminiscent of a TSWV infection, was received by the Netherlands Inspection Service for Horticulture. Preliminary serological tests with available antisera to a wide range of different orthotospoviruses were negative and indicated the possible emergence of a new orthotospovirus species. Here data is presented on the biological, serological and molecular characterization of this Als-2000 virus isolate to support the proposal of a new tentative orthotospovirus species for which the name alstroemeria yellow spot virus (acronym AYSV) is proposed. The nucleotide sequence of the entire genome was elucidated using a simultaneous conventional and NGS cloning/ sequencing approach.

Ribonucleocapsid purification, SDS-PAGE and serological analysis
Ribonucleocapsid protein (RNP) was purified from Als-2000-infected N. benthamiana plants [47] and subjected to 25-45% CsSO 4 gradient centrifugation. RNP fractions collected were used to immunize rabbits for the production of a polyclonal antiserum. Prior to this, RNP fractions were subjected to a quality/purity check, in comparison to INSV and TSWV RNP preparations, on a 12% SDS-PAGE after denaturation in Laemmli buffer [48]. Rabbits were immunized 2-3 times with 50-100 μg of purified RNP emulsified with incomplete Freund´s adjuvant (1:1, w/v), during intervals of two weeks. Serum was collected two weeks after the last injection and used for the purification of IgG and an the preparation of an IgG-alkaline-phosphatase conjugate. For orthotospovirus differentiation double-antibody sandwich-enzyme-linked immunosorbent assays (DAS-ELISA) were performed using a panel of polyclonal antisera raised against Als-2000, TSWV, TCSV, GRSV, IYSV and TYRV [45][46][47]. Sap of infected and healthy N. benthamiana plants was diluted 1:30 and used as the antigen source and negative control, respectively.

RNA purification and RT-PCR
Viral RNA was purified either from purified RNPs using 1% SDS and hot phenol extraction followed by ethanol precipitation [9] or from infected plant material using Trizol buffer according to the manufacturer's instructions (Invitrogen).
In a first approach to obtain Als-2000 viral sequences only a single primer "Asian Termini" (AT; 5΄-CCC GGA TCC AGA GCA ATC GAG -3΄, containing a BamHI restriction site (underlined) and matching the 3' terminal sequences of Asian orthotospoviruses (italics), was used for firststrand cDNA synthesis using M-MLV Reverse Transcriptase (Promega), and a subsequent PCR amplification [46] using Phusion High-Fidelity DNA Polymerase (New England BioLabs). Since the sequence of the amplicon obtained using the AT primer showed highest homology with the NS M gene of IYSV, a second approach was applied in which Als-2000 sequences were additionally RT-PCR amplified using primers IY1 (5′-CCCGA GGA TCC ATG GCT ACC GTT AGGG -3′) and IY2 (5′-CCCGA GGA TCC AAA TTA ATT ATA TCT ATC TTT CTT GG-3′) corresponding to IYSV N gene sequences (in italics) combined with the universal hairpin primer (UHP) (5'-CACT GGA TCC TTT TGT TTT TGT TTT TTG -3') as described before [45,46]. Based on newly obtained sequences for Als-2000, specific primers were designed to verify/ obtain additional sequences to complete the entire S RNA nucleotide sequence. Primers Als2000-NS M -R (5'-GCT TGG TTT TCT TTC TTT TTC CTT C-3') and Als2000-G-R (5'-CCT CCT AAA ACA TAT GAC TTTCC-3') were used to amplify the 5' ends of the NS M respectively G N /G C glycoprotein precursor ORFs in combination with the AT termini primer. In analogy to the approach for the S RNA segment, new and specific primers were designed to verify and/or obtain sequences for the intergenic region of the M RNA segment. All PCR fragments obtained were cloned into pJET blunt1.2 vector (Fermentas) and their sequences determined (Eurofins).

Sequence analysis of S and M RNA segments
Contigs for the M and S RNA segments were assembled from obtained sequences using DNAStar software and used to search the GenBank database using BLASTn and BLASTp. Multiple sequence alignment and phylogenetic analysis were performed using MEGA 7 software [49]. The G N -G C sequence was analyzed for the presence of N-and O-linked glycosylation sites by prediction using the NetNGlyc 1.0 and NetOGlyc 3.1 Servers (http://www.cbs.dtu.dk/ servi ces/NetNG lyc/; http://www.cbs.dtu.dk/servi ces/NetOG lyc-3.1/). Potential signal cleavage sites within the glycoprotein were predicted using the SignalP 4.1 server (http:// www.cbs.dtu.dk/servi ces/Signa lP/). Transmembrane regions within the G N -G C precursor protein were predicted using the TMHMM Server v.2.0 program (http://www.cbs.dtu.dk/servi ces/TMHMM /).

Sequences analysis of Als-2000 by Next Generation Sequencing
For the sequencing of the entire Als-2000 genome using next generation sequencing (NGS) technology, total RNA was isolated from systemically infected N. benthamiana leaves with the Qiagen RNeasy plant mini kit according the manufacturers' procedures (Qiagen). Samples were prepared for HighSeq sequencing according to Verbeek et al. [50]. Sequence data were analyzed using the CLC Genomics Workbench 8.5.1 (Qiagen) as earlier described [51]. To obtain and verify the full length sequence of Als-2000 the 5'-and 3'-ends of the viral RNA segments and regions from which only few reads were obtained, were RT-PCR amplified with specific primers and directly sequenced by Sanger sequencing (Macrogen).

Thrips transmission
To determine the potential of Thrips tabaci for transmission of Als-2000, 10 first and early second instar larvae of T. tabaci were given an acquisition period of 24h on Als-2000 infected Emilia sonchifolia and next transferred to healthy E. sonchifolia plants. Plants were kept under cages at room temperature and monitored for symptom expression for two weeks. The presence of Als-2000 in plants was tested by RT-PCR using primers IY1 and IY2. Fragments of expected sizes were collected and the AYSV sequence confirmed by sequence analysis.

Host range analysis
From 31 plant species mechanically inoculated with Als-2000, from here onwards referred to as alstroemeria yellow spot virus (AYSV), 17 species did not show any sign of infection, 9 species only showed local symptoms while 5 species displayed local and systemic symptoms ( Table 1). Back-inoculation of the isolate from N. benthamiana on the alstroemeria cultivar Dimention did not induce any symptoms. The data indicated that the virus had a relatively narrow host range and among the plants tested only dicots showed local or systemic symptoms.

N protein electrophoresis and DAS-ELISA
To comparatively analyze AYSV and serologically differentiate this virus from other orthotospoviruses, AYSV RNPs were purified and analyzed by SDS-PAGE. A major band of approximate M W 30 kDa migrated at a position in between the INSV and TSWV N proteins and likely represented the AYSV N protein (Figure 2, left panel). Polyclonal antibodies produced against the AYSV RNPs were used to analyze its serological relationship to a set of orthotospoviruses. In DAS-ELISA, only leaf extracts of AYSV-infected plants reacted positively with polyclonal antiserum against AYSV, while no positive reaction was observed with antisera raised against other orthotospovirus species (Figure 2, right panel).

Sequence analysis and characterization of the S RNA, M RNA and L RNA
The entire sequence of AYSV S RNA consisted of 2734 nucleotides (nt) with 5′ terminal (72 nt) and 3′ terminal (71 nt) untranslated regions (UTR) that contain the highly conserved 8 terminal nucleotides present in all orthotospovirus RNA genome segments. AYSV terminal sequences show almost entire base pair complementarity of the first 18 nt with only 3 base pair mismatches (5'-AGA GCA AUCGNNGNAUAA-3'). This complementarity continues, although to a lesser extent, to ~ 60 nucleotides of the UTR sequences and assists in the panhandle formation of the S RNA segment.
The TYRV S RNA contains a vORF starting with an AUG at nucleotide position 73 and terminating with a UAG codon at position 1404, coding for 443 residues of the NS S protein with a predicted molecular mass of 50.1 kDa. The Upon multiple sequence alignment of the AYSV N protein to other orthotospoviral N proteins, AYSV N showed highest identity (69.5 %) with polygonum ringspot virus (PolRSV). This homology is well below the threshold for demarcation of new orthotospovirus species and indicated that AYSV tentatively presents a new distinct species. In a phylogenetic analysis based on the N protein, AYSV clusters together with members initially found in Europe and Asia that belong to the so-called Eurasian clade. A similar phylogenetic analysis using AYSV NS S revealed 88.3 % identity to that of TYRV, and supports the relationship of AYSV to the Eurasian clade ( Figure 3).
The complete nucleotide sequence of AYSV M RNA consists of 4797 nt, and contains a 5'-(62 nt) and 3'-UTR (51 nt) that is able to form a panhandle structure due to sequence complementarity. The UTRs show almost entire sequence complementarity for 29 nt with only 3 mismatches (5'-AGA GCA AUC GGU GCAANAAU CAA AUNUNU-3'). The vORF starts with an AUG at nucleotide position 63 and terminates with an UAA codon at position 989, coding for the NS M protein with a predicted molecular mass of 34.2 kDa. The AYSV NS M protein shows conserved motifs that belong to the 30K movement protein superfamily [52,53], including D 155 , G 211 and P/D-L-X (amino acid residues D 102 S 103 L 104 ), and motifs for a phospholipase A2 catalytic site and PLA2 (amino acid residues CMQLNLTS [197][198][199][200][201][202][203][204] ).
The vcORF codes for the glycoprotein precursor and starts with an AUG codon at nucleotide position 4746 and terminates with UGA at nucleotide position 1300. The M RNA encoding v-and vc-ORFs flank an IGR of 310 nt that runs from nucleotide position 990 to 1299 and, like the S RNA, consists of nucleotide stretches highly rich in A-and U residues. Folding predictions of the IGR sequence reveals the formation of a stable hairpin structure (data not shown).
Using BLASTp, the highest identity of AYSV NS M was found with TYRV NS M (88.7%) whereas its glycoprotein precursor (GP) showed highest identity with that of TYRV GP (81.0%). Phylogenetic analysis based on the NS M protein and the glycoprotein precursor support a classification of AYSV within the Eurasian clade ( Figure 3).
The surface glycoprotein of AYSV was analyzed and compared to those of other Eurasian orthotospovirus species for the presence of N-, O-Glycosylation, peptide cleavage sites and transmembrane domains (Figure 4). Two putative N-glycosylation sites, Asn 319 and Asn 575 and eight O-glycosylation sites i.e. Thr 80 , Thr 82 , Thr 84 , Thr 99 , Thr 100 , Thr 101 , Thr 102 and Thr 105 were predicted. Two signal peptide sequences with predicted cleavage sites were found at amino acid positions Leu 23 and Leu 462 . Five (hydrophobic) transmembrane domains were found in the G N -G C precursor. Both of the predicted cleavage sites located within the first and fourth transmembrane domains, whereas none of the predicted glycosylation sites mapped within transmembrane domains (Figure 4). A comparison to the G N -G C precursor Based on the predicted cleavage sites, the AYSV G N -G C precursor is processed into the mature G N and G C glycoproteins of 50.3 and 77.4 kDa, respectively (Figure 4). The AYSV G N -G C precursor does not contain a Arg-Gly-Asp (RGD) sequence, a motif earlier found in the TSWV precursor and putatively involved in virus cell attachment [9].
The entire sequence of AYSV L RNA consists of 8865 nt with a 5′ terminal (214 nt) and 3′ terminal (32 nt) UTR that contain the highly conserved 8 terminal nucleotides present in all orthotospoviral RNA segments. AYSV terminal sequences show almost entire base pair complementarity for the first 16 nt with only one base pair mismatch (5'-AGA GCA AUCNAGC AAC -3').
The L RNA contains a single ORF of 8619 nt in the vcstrand and coding for the putative viral RNA-dependent RNA polymerase (RdRp, L protein) of 2873 aa with a predicted molecular mass of 316.0 kDa. The vcORF starts with an AUG at nucleotide position 8833 and terminates with an UAA stop codon at nucleotide position 212. The L gene of AYSV shares highest aa (87.3 %) identity with the one from TYRV. Analysis of the AYSV L protein revealed the presence of all conserved motifs earlier reported and indicative for the RdRp [54,55]. The highly conserved region, shared among all characterized orthotospoviral L proteins, is found at aa position 1,275-1,559 of the AYSV L protein. This region contains the essential motifs of viral RdRp proteins from negative strand RNA viruses and includes motif A (aa 1,365-1,378), motif B (aa 1,452-1,471), motif C (aa 1,490-1,499), motif D (aa 1,534-1,544), motif E (aa 1,547-1,559), and motif F (aa 1,275-1,301) [54,55]. The AYSV RdRp motif sequences are 100% identical to those of TYRV. Phylogenetic analysis of the RdRp shows a clustering of AYSV with TYRV, PolRSV, hippeastrum chlorotic ringspot virus (HCRV) and IYSV (Figure 3).
A summary and comparative analysis of the AYSV L, M and S RNA sequence to those of other Eurasian species is shown in Table 2. While the L RNA sequence has not been elucidated for all Eurasian species, the structural features of the M and S RNA of AYSV show highest similarity to those from PolRSV, the closest related species, with the exception of a differently sized IGR within the S segment (AYSV-IGR 478 vs. PolRSV-IGR 183 ).
The complete L, M and S RNA sequences of AYSV are available under GenBank accession no. MF469033, MF469034 and MF469035 respectively

Vector transmission
Healthy E. sonchifolia plants exposed to T. tabaci, that in their larval stages were given access to AYSV-infected plants, revealed chlorotic spots similar to those observed earlier during the host range analyses ( Table 1). The presence of AYSV was confirmed and successful transmission by viruliferous thrips was demonstrated by RT-PCR amplification of (partial) AYSV N gene sequences from infected leaf material and subsequent sequence analysis.

Discussion
Here a tentative new orthotospovirus species has been identified and characterized from alstroemeria, and for which the name alstroemeria yellow spot virus (AYSV) has been proposed. The virus causes yellow spots on alstroemeria, has a narrow experimental host range limited to a few dicots and is transmitted by the vector T. tabaci. A polyclonal antiserum has been raised against AYSV RNPs and only reacts positively with the homologous virus and not with any other established orthotospovirus species tested. The nucleotide sequence of the genomic L, M and S RNA segments has been elucidated (8619, 4797 and 2734 nt, respectively) and contains 5 ORFs coding for six mature proteins i.e., the L protein (viral RdRp, 316 kDa), the cell-to-cell movement protein (NS M , 34.2 kDa), the precursor (127,7 kDa) to the glycoproteins G N and G C , the suppressor of RNA silencing (NS S , 50.1 kDa) and the nucleocapsid protein (N, 30.1 kDa). The AYSV sequence shows all structural features typical for orthotospovirus RNA genomes with highest homology between its proteins and those from TYRV. Phylogenetic analysis indicates that AYSV, together with PolRSV, TYRV, HCRV and IYSV, belongs to a growing clade of ◂ orthotospoviruses for which the geographical origin seems to be based in the Middle-East [46]. Although the origin of the virus-infected alstroemeria plant could not be indicated with certainty, it was provided by a grower that at that time had simultaneously received shipments of alstroemeria plants from Colombia and Iran. Considering the relation of AYSV to the orthotospoviruses from the Middle-East cluster makes it most likely that Iran is the site of origin of the collected sample.
With the identification of this new orthotospovirus species, the family of Tospoviridae is still expanding and has currently reached the number of 30 tentative and established species ( [2] and, this report). While this new species has only a small (experimental) host range, its natural host species so far is limited to alstroemeria. With the availability of a specific polyclonal antiserum to AYSV, field surveys could help to identify other natural (cultivated and non-cultivated) plant hosts of this virus. Considering that nothing is known on the prevalence and distribution of AYSV in agricultural and horticultural crops and ornamentals these surveys will also indicate the economic importance of this virus.
During this study a combination of conventional and NGS sequencing approaches has helped to elucidate the entire AYSV genome sequence. The speed and complete data sets Fig. 4 Topology of the glycoprotein precursor of AYSV in comparison to the ones from other Euroasian orthotospoviruses. Predicted glycosylation sites are indicated with arrows op top. Hydrophobic domains representing signal sequences and transmembrane domains are indicated in grey, and their locations within the precursor with the position of amino acids below. Predicted cleavage sites for maturation of the precursor into the mature G N and G C glycoproteins are indicated with scissors with the amino acids sequence of the cleavage sites above it (TEV-L and SLA-L, respectively) Table 2 Comparison of structural features of the L, M and S RNA segments of AYSV to those of other Eurasian orthotospoviruses obtained by the NGS approach, combined with lowering costs, favor usage of NGS for the elucidation of viral RNA/ DNA genomes. However, during this study and our earlier sequence analysis of chrysantemum stem necrosis virus (CSNV) [51] we have observed that the NGS approach may not always provide the entire 5' and 3' -terminal sequences and sufficient reads to span the entire A-and U-rich intergenic region. Assembly of the latter reads into one contig is also hampered by the fact that they are relatively short and subsequently fail to properly assemble due to sequence stretches rich in A-and U-residues. Verification of the termini and intergenic region thus still needs to be done by conventional RT-PCR / 5' and 3' RACE amplification and Sanger sequence analysis of the amplicons.


extfieldsofstudy:
Biology
Computer Science

text:
Identifying protein subsets and features responsible for improved drug repurposing accuracies using the CANDO platform

Drug repurposing is a valuable tool for combating the slowing rates of novel therapeutic discovery. The Computational Analysis of Novel Drug Opportunities (CANDO) platform performs shotgun repurposing of 3,733 drugs/compounds that map to 2,030 indications/diseases by predicting their interactions with 46,784 protein structures and relating them via proteomic interaction signatures. The accuracy of the CANDO platform is evaluated using our benchmarking protocol that assesses indication accuracies based on whether or not pairs of drugs associated with the same indication can be captured within a certain cutoff, which is a measure of the drug repurposing recovery rate. To identify subsets of proteins that exhibit the same therapeutic effectiveness as the full set, groups of 8 proteins were randomly selected and subsequently benchmarked 50 times. The resulting protein sets were ranked according to average indication accuracy, pairwise accuracy, and coverage (count of indications with non-zero accuracy). The best 50 subsets of 8 according to each metric were progressively combined into supersets after each iteration and benchmarked. These supersets yield up to 14% improvement in benchmarking accuracy, and represent a 100-1,000 fold reduction in the number of proteins relative to the full set. Protein supersets optimized using independent compound libraries derived from the full library were cross-tested and were shown to reproduce the performance relative to using all 46,784 proteins, indicating that these reduced size supersets are broadly applicable for characterizing drug behavior. Further analysis revealed that sets comprised of proteins with more equitably diverse ligand interactions are important for describing drug behavior. Our work elucidates the role of particular protein subsets and corresponding ligand interactions that play a role in computational drug repurposing, and paves the way for the use of machine learning approaches to further improve the accuracy of the CANDO platform and its repurposing potential. Author summary Drug repurposing is a valuable approach for ameliorating the current problems plaguing drug discovery. We introduce a novel protein subset analysis pipeline that allows us to elucidate features important for drug repurposing accuracies using the Computational Analysis of Novel Drug Opportunities (CANDO) platform. Our platform relates drugs based on the similarity of their interactions with a diverse library of proteins. We subjected all proteins in the platform to a splitting and ranking protocol that ranked protein subsets based on their benchmarking performance. Further analysis of the best performing protein subsets revealed that the most useful proteins for describing how small molecule compounds behave in biological systems are those that are predicted to interact with a structurally diverse range of ligands. We hypothesize that this is a consequence of the multitarget nature of drugs and, conversely, the implied promiscuity of proteins in biological systems. These results may be used to make drug discovery more accurate and efficient by alleviating some of its bottlenecks, bringing us one step further in better understanding how drugs behave in the context of their environments.

Common strategies in drug discovery include forward pharmacology [1] and rational 2 drug design [2]. In the former, a library of compounds is screened, typically in high 3 throughout manner, for certain phenotypic effects in vitro. In the latter, compounds are 4 virtually screened against a predetermined biological target and high confidence hits are 5 then assayed for a desired modulation. In both cases, the hits obtained are then 6 assessed for effectiveness in vivo and proceed to clinical trials for eventual FDA 7 approval if successful at each stage. This iterative process can cost billions of dollars 8 and up to 15 years per drug [3]. These approaches do not consider the promiscuity of 9 approved drugs in the context of indications/diseases within living systems (evidenced 10 by side effects present for all small molecule therapies [4,5]), dooming many novel 11 therapeutics to fail. With the second-leading cause of putative drug attrition being 12 adverse reactions [6], there is great utility in finding new uses for already approved 13 drugs, which is formally known as drug repurposing or repositioning [7,8]. 14 We have developed the Computational Analysis of Novel Drug Opportunities 15 (CANDO) platform [9][10][11] to address these drug discovery challenges. One fundamental 16 tenet of CANDO is that drugs interact with many different proteins and pathways to 17 rectify disease states, and this promiscuous nature is exploited to relate drugs based on 18 their proteomic signatures [9,[12][13][14][15]. These signatures are typically determined via 19 virtual molecular docking simulations that are applied to predict compound-protein 20 interactions on a proteomic scale. Using a knowledge base of known drug-indication 21 approvals/associations, we can identify putative drug repurposing candidates for a 22 particular indication based on the similarity of their proteomic interaction signatures to 23 all other drugs approved for (or associated with) that indication. When a particular 24 indication does not have any approved drug, the library of human use compounds 25 present in CANDO is screened against the tertiary structures of all relevant and 26 tractable proteins obtained by x-ray diffraction or homology modeling from a particular 27 organismal proteome to suggest new treatments that maximize binding to the 28 disease-causing proteins and minimize off-target effects. High-confidence putative drug 29 candidates generated by CANDO using both approaches have been prospectively 30 validated preclinically for a variety of indications, including dengue, dental caries, 31 diabetes, hepatitis B, herpes, lupus, malaria, and tuberculosis, with 58/163 candidates 32 yielding comparable or better therapeutic activity than standard 33 treatments [12,13,16,17]. 34 To date, putative drug candidates generated by CANDO have been based on simple 35 comparison metrics, primarily the root mean square deviation (RMSD) of the binding 36 scores present in a pair of drug-proteome interaction signatures. Our platform is 37 evaluated using a benchmarking method that assesses per indication accuracies based 38 on whether or not other drugs associated with the same indication can be captured 39 within a certain cutoff in terms of similarity to a particular drug approved for that 40 indication. Incorporating machine learning, which is continuing to prove its utility in 41 many aspects of biomedicine [18][19][20] including drug discovery and repurposing [21,22], 42 August 24, 2018 2/17 into the CANDO platform to increase benchmarking accuracies and therefore its 43 predictive power is of importance. Various algorithms can be incorporated (for example, 44 neural networks, support vector machines, and decision trees), but the well documented 45 issues described by the curse of dimensionality [23,24]  OpenBabel FP4 fingerprinting method [37], resulting in a structural similarity score.

103
The score that populates each cell in the compound-protein interaction matrix is the 104 maximum value of all of the possible binding site scores times the structural similarity 105 scores of the associated ligand and the compound.

106
Benchmarking protocol and evaluation metrics

107
The compound-compound similarity matrix is generated using the root mean square 108 deviation (RMSD) calculated between every pair of compound interaction signatures 109 (the vector of 46,784 real-value interaction scores between a given compound and every 110 protein in the library). Two compounds with a low RMSD value are hypothesized to 111 have similar behavior [9-11, 13, 15]. For each of the 1,439 indications with two or more 112 associated drugs, the leave-one-out benchmark assesses accuracies based on whether 113 another drug associated with the same indication can be captured within a certain 114 cutoff of the ranked compound similarity list of the left-out drug. This study primarily 115 focused on a cutoff of the ten most similar compounds ("top10"), the most stringent 116 cutoff used in previous publications [9-11, 13, 15]. The benchmarking protocol calculates 117 three metrics to evaluate performance: average indication accuracy, (top10 in this study) and d is the number of drugs approved for that given indication.

122
Pairwise accuracy is the weighted average of the per indication accuracies based on how 123 many drugs are approved for a given indication. Coverage is the count of the number of 124 indications with non-zero accuracies within the top10 cutoff.

125
Superset creation and benchmarking 126 The 46,784 proteins in the CANDO platform were randomly split into 5,848 subsets of 8 127 and subsequently benchmarked using the method described above. The size of 8 was 128 selected because offered the widest range of benchmarking values (relative to larger 129 sizes), reduced the computational cost of the experiments (relative to smaller sizes that 130 increase the number of individual benchmarks that need to be evaluated), divided into 131 46,784 evenly, and also provided adequate signal for the multitargeting approach to work 132 according to our prior studies [12]. A total of 50 iterations were performed that resulted 133 in 292,400 benchmarking experiments. Each subset was then ranked according to top10 134 average indication accuracy, pairwise accuracy, and coverage. The fifty best performing 135 subsets from each ranking criterion (average indication accuracy, pairwise accuracy, and 136 coverage) were progressively combined into supersets and benchmarked after each of the 137 50 iterations of the splitting and ranking protocol. The subsets were nonredundantly combined such that if a given protein was represented in the best performing subsets The best performing protein subsets and supersets were further analyzed to elucidate 160 the protein feature(s) responsible for increased benchmarking performance. The protein 161 subsets and supersets were analyzed based on four criteria: organismal source, structure 162 source (x-ray diffraction or modeling), fold space (based on the CATH classification of 163 proteins [38]), and interacting ligand structure distributions. The subsets and supersets 164 were analyzed by counting the specific organisms to which the proteins belonged to see 165 if any were over or underrepresented in the best and worst performing sets. Similarly, 166 the subsets and supersets were analyzed to see if structures obtained via a specific 167 source, x-ray diffraction or modeling, were differentially represented in the best and 168 worst performing sets. Fold assignments were made to each protein in the subsets and 169 supersets which were again analyzed for differential representation of specific protein 170 folds. Finally, since our compound-protein interaction scoring method utilized the 171 structural similarity of each drug/compound to the ligand co-crystallized with the 172 protein (see previous section), we analyzed these ligands for differential representation. 173 Each co-crystallized ligand in the COFACTOR database of template binding sites were 174 clustered at various distances (0.1 to 0.9 with increments of 0. performance for each superset within a given iteration tends to gradually increase as the 212 number of iterations increases (Fig 1). Ranking subsets using coverage is overall the  Sorting the supersets by size reveals that at least 80-120 proteins are required to 223 reach optimal benchmarking performance (Fig 2). Nonredundantly combining the worst 224 performing subsets of 8 into subpar sets demonstrates worse performance than the 225 control value for the average indication accuracy based on using the full library, with 226 the mean values of these subpar set distributions being below the acceptable 5% 227 threshold of 11.6% for all sizes benchmarked. The random set and subpar set 228 distributions begin to converge toward the average indication accuracy control value as 229 size increases (Fig 2).  being the best ranking criterion through ten iterations, after which coverage is superior, especially when measuring the coverage metric. This result demonstrates that the splitting and ranking protocol can produce supersets with benchmarking performance superior to using the full protein library by combining the best performing subsets with a vastly reduced number of proteins (100 to 1,000 fold reduction in size), further suggesting that specific groups of proteins are relatively more useful for drug repurposing accuracy.

Cross-testing with independent compound libraries 234
For the independent compound library experiment, average indication accuracy was 235 chosen as the ranking criterion because it performed the best in all three metrics 236 through ten iterations in the superset experiment (Fig 1). Box and whisker plots for 237 each compound library show the spread of the benchmarking performance for each 238 metric generated from the protein supersets obtained through the splitting and ranking 239 protocol of their complementary compound library (Fig 3). Supersets less than size 80 240 were excluded because there is a minimum number of protein features required to reach 241 optimal benchmarking performance (Fig 2). The control value for the given metric fell 242 within the inter-quartile range or below in 26 cases (87%), while it was within the upper 243 quartile in the remaining four cases (13%), indicating that these values always fell within 244 the range of the accuracy/coverage distributions. The control value was an extreme Superset, subpar set, and random set average indication accuracies sorted by size. Average indication accuracies are shown for the supersets (blue) generated using the best subsets ranked by the same metric. The line traces the mean score for each size with the bars indicating one standard deviation for the distribution. Subpar sets (orange) are the combinations of the worst performing subsets ranked by average indication accuracy. Randomly selected protein sets (green) of each size were also generated and benchmarked. The control value based on using the full protein library (dashed black at 12.3%) and an acceptable 5% threshold (dotted black at 11.6%) are plotted for reference (i.e., any protein set that benchmarks within 95% of the control value is considered acceptable). For the random sets and supersets, the performance in terms of average indication accuracy begins to plateau around 80-120 proteins. The supersets begin to slightly decline in performance after 32-33 subsets (256-264 proteins). The mean subpar set accuracies at each size all fall outside the 5% acceptable threshold, while the superset distributions are well above the control value with as few as five subsets. The difference between the superset and subpar set performance suggests that there is a particular distribution of features within their proteins that is correlated with benchmarking performance. supersets to describe compound behavior more effectively than the full protein library. 247

Ligand clustering and feature-based creation of protein libraries 248
The protein subsets and supersets were analyzed based on four criteria to elucidate the 249 feature(s) responsible for benchmarking performance: organismal source, structure 250 source (x-ray diffraction or modeling), fold space, and interacting ligand structure 251 distributions. There were no significant correlations found for the first three criteria; no 252 organism(s) or fold(s) was consistently represented in the best performing sets, while  The solid blue and solid orange points are averages of 50 best and worst subsets, respectively. Dashed lines represent an example of a superset (blue) and a subpar set (orange) which are nonredundant combinations of the best and worst performing subsets respectively. The control set (dashed black), representing the full protein structure library, falls in between the superset and subpar set. The black diamonds indicate that the distribution of counts at that cluster rank between the best and worst performing subsets, assessed using Welch's T-test, is significant (p-value < 0.05). The subsets and supersets with the best performance demonstrate a more equitable distribution of interactions among ligand clusters as opposed to the worst performing subsets and subpar sets, indicating that using multitargeting proteins to compose our structure libraries yields superior benchmarking performance.
Based on the subset analysis, it was hypothesized that proteins having a more 276 diverse and equitable distribution of interactions will benchmark better. Protein protein, we applied a upper cutoff filter for the number of ligand clusters allowed (Fig 284  5); all libraries were limited to size 80, which is the minimum required to reach optimal 285 performance based on data in Fig 2. There is an optimal upper limit of ≈20-40 ligand 286 clusters with regard to the benchmarking, indicating that using too many ligand 287 clusters to describe a protein is undesirable for characterizing drug/compound behavior. 288 Based on the data in Fig 5, we created protein libraries of various sizes using a upper 289 cutoff of 40 (Fig 6). We begin to consistently recapture the benchmarking performance 290 of the full library (within 5% error) with as few as 60 proteins. Libraries composed of The splitting and ranking protocol was originally intended to find a protein subset that 300 benchmarked as least as well as the full set. The improvement of the benchmarking 301 performance is an encouraging sign for incorporating machine learning in the CANDO 302 platform in the future, and discovering how more complex weighting and relating of 303 proteins contribute to drug repurposing accuracy, which is difficult to do with simple 304 RMSD calculations. The smaller-sized protein libraries generated as part of this study, 305 representing a 100 to 1,000 fold reduction in size, will be more conducive to machine 306 learning. Feature reduction through the use of neural network based auto-encoders or 307 principal component analysis will provide an important contrast to our proposed 308 method.

309
The independent compound library experiment demonstrated that optimized protein 310 sets based on a particular library were capable of therapeutically characterizing a 311 completely different one, indicating that these supersets are generalizable. In other 312 words, if a new drug/compound is added to the CANDO putative drug library, these 313 reduced size supersets are likely able to describe its behavior at least as well as using 314 every protein available. In addition to facilitating machine learning, our findings suggest 315 a greatly reduced time required to generate new proteomic interaction vectors, which is 316 particularly important if the program/protocol of choice for generating interactions is 317 computationally expensive. Any repurposing candidates suggested from using the 318 supersets are on average more clinically relevant as they were able to recapture drug 319 behavior more accurately than using the full protein library in a statistically significant 320 manner.

321
The knowledge that ligand cluster interaction diversity is the key requisite feature in 322 describing drug behavior may also lead to optimization strategies other than the 323 computationally intensive splitting and ranking protocol used in this study. The . A minimum of two ligand clusters was required to avoid the trivial case of only one cluster mapped (with a variance of zero). A upper ligand cluster count cutoff was applied to these protein libraries to determine the effect on benchmarking performance. The blue line traces the average indication accuracy distribution mean from benchmarking 50 libraries at each cutoff using the top ranked proteins. Similarly, the orange line traces the average indication accuracy distribution mean from benchmarking 50 libraries at each cutoff using the highest variance proteins. The bars indicate one standard deviation for the distribution. Using too small (< 20) or too large (> 40) of a cutoff results in suboptimal benchmarking performance. The high variance libraries result in average performance far below the acceptable 5% error range, with the cluster cutoff seemingly having no effect on average indication accuracy as all cutoffs produced comparable distributions. All average indication accuracies produced using a upper cutoff of 20-40 ligand clusters were within the acceptable 5% error range, with the upper cutoff of 40 ligand clusters producing the greatest mean accuracy. This result indicates that there is an optimal range of ligand cluster interactions to best describe therapeutic behavior. therapeutically relevant if the proteins used to describe the behavior of a compound are 333 diverse in terms of the structures of the ligands which interact with their binding sites. 334 Protein libraries with fewer predicted ligand cluster binding partners yield much worse 335 performance than those consisting of proteins interacting with a more structurally 336 diverse range of ligands. Coupling this with the finding that there is a minimum 337 number of proteins required to reach optimal benchmarking accuracies (Fig 2), which 338 was also observed by us previously [10,11], drugs should realistically be described in the 339 Using too small of a size (< 60) results in suboptimal benchmarking performance. Creating libraries from proteins with the highest variance results in performance on average far below the acceptable 5% error range, although size does have a positive correlation with performance with these high variance sets. This result reiterates that there is a minimum number of proteins required to reach optimal benchmarking performance and that proteins with high variance in their ligand cluster signatures are far superior for describing drug/compound behavior.
context of their multitarget nature, treating both small molecule compounds and 340 proteins promiscuously, as in biological systems [43][44][45]. However, using libraries of 341 proteins with too many diverse interactions in the CANDO platform also leads to 342 suboptimal performance. We hypothesize this can be attributed to two factors: 1) 343 spreading a compound interaction signature across too many (50 or more) ligand 344 clusters can potentially dilute the therapeutic signal relative to the promiscuity of the 345 corresponding proteins; or 2) these proteins are not therapeutically relevant and are 346 therefore not useful for specifically describing drug behavior.   Visualization of the best and worst protein types for benchmarking performance based on ligand cluster signatures. Centroids of the top five ligand clusters from the signatures of each protein are depicted. The percent of interactions belonging to each ligand cluster are next to their respective centroid. Surface representations of the proteins were made using Chimera [42] with the predicted binding site residues from COFACTOR for each ligand shown (excluding the smaller ligands in C) colored in blue. A) Alkaline serine protease KP-43 from Bacillus subtillus: the top five ligand clusters account for 46.3% of the total interactions with the distribution between them being relatively equitable. B) SET domain of human histone-lysine N-methyltransferase: only two ligand clusters are predicted to interact, with one having over 98% of the total interactions. C) Human STE20-related kinase adapter protein beta: the ligand cluster signature is too promiscuous with the top five ligand clusters accounting for only 18.3% of the total interactions; the remaining sixteen ligands surrounding the protein account for 28.7%, which combined with the top five clusters is as much as the 46.3% of the total interactions shown in A from only five clusters. Subsets and supersets composed of proteins similar to A outperform those composed of proteins similar to B and C in benchmarking, indicating that moderately promiscuous proteins with equitable ligand cluster signatures are the best therapeutic descriptors.

348
We have developed an integrated pipeline that allows for the elucidation of proteins and 349 their features that are important for benchmarking in the CANDO platform, and 350 therefore important for drug repurposing. We are able to reproduce the performance of 351 the complete CANDO protein structure library with orders of magnitude fewer proteins, 352 allowing for more rapid candidate generation when evaluating new putative drug 353 libraries or any other changes to the platform. We discovered that moderately


extfieldsofstudy: Biology Medicine

text:
Differential gene transcription across the life cycle in Daphnia magna using a new all genome custom-made microarray

Background Unravelling the link between genes and environment across the life cycle is a challenging goal that requires model organisms with well-characterized life-cycles, ecological interactions in nature, tractability in the laboratory, and available genomic tools. Very few well-studied invertebrate model species meet these requirements, being the waterflea Daphnia magna one of them. Here we report a full genome transcription profiling of D. magna during its life-cycle. The study was performed using a new microarray platform designed from the complete set of gene models representing the whole transcribed genome of D. magna. Results Up to 93% of the existing 41,317 D. magna gene models showed differential transcription patterns across the developmental stages of D. magna, 59% of which were functionally annotated. Embryos showed the highest number of unique transcribed genes, mainly related to DNA, RNA, and ribosome biogenesis, likely related to cellular proliferation and morphogenesis of the several body organs. Adult females showed an enrichment of transcripts for genes involved in reproductive processes. These female-specific transcripts were essentially absent in males, whose transcriptome was enriched in specific genes of male sexual differentiation genes, like doublesex. Conclusion Our results define major characteristics of transcriptional programs involved in the life-cycle, differentiate males and females, and show that large scale gene-transcription data collected in whole animals can be used to identify genes involved in specific biological and biochemical processes. Electronic supplementary material The online version of this article (10.1186/s12864-018-4725-7) contains supplementary material, which is available to authorized users.

Background
Unravelling the link between genes and environment across a life-cycle is a challenging goal. The full exploration of this link requires model organisms with wellcharacterized ecological interactions in nature, tractability in the laboratory and available genomic tools [1]. In Drosophila, the vast majority of genes are transcribed during the first embryonic stages due to the contribution of both maternally-inherited and zygote-produced mRNAs, while adult females and males display specific transcribed genes, consistent with the known transcriptional profiles of ovaries and testis [2,3]. Embryogenesis in the worm Caenorhabditis elegans elicits major transcriptional changes, particularly during the early and late stages, with a significant lower number of gene expression changes in mid-embryogenesis [4]. Males of C. elegans show a large number of highly expressed, malespecific genes, most of which being undetectable in any hermaphrodite stage [4]. Similar transcription gene patterns during embryogenesis and adult stages have been reported for many other species of sponges and insects [5,6]. In the above mentioned species, the observed high transcription of genes during embryogenesis are associated with morphogenesis, differentiation, and development, and larval or adult stages show different clusters of gene transcription often related with particular tissues and organs.
The waterflea Daphnia magna Strauss is a very useful species for both environmental and basic biology studies [7,8]. D. magna can be found in freshwater ecosystems all around the world, playing a relevant ecological role in the food web, both as a primary consumer and as a food source for predators, both vertebrates and invertebrates [9]. D. magna can be maintained in the lab as genetically uniform clonal lines through asexual parthenogenic reproduction, hereby reducing genetic variability and enabling the parallel analysis of functional and fitness changes in the same genotype in multiple environmental conditions [10]. D. magna and its close relative D. pulex are widely used as indicators of environmental quality [11], as models in evolutionary biology, and in the study of adaptive responses to environmental changes [12][13][14][15][16].
Daphnids genomes have recently received much attention from the scientific community, with the publication of two full genomes (D. pulex and D. magna) and two more species having their genomes partially sequenced [1,[17][18][19][20]. The specific function of many genes, however, even in well studied organisms, is incomplete. It is not unusual to find a third of the genes in any genome to lack functional annotation [21,22]. This lack of functional annotation is higher in newly sequenced organisms that often contain lineage-specific genes lacking homology with other assembled genomes [22,23]. This is the case for D. magna and its close relative D. pulex, [1,20,23,24]. In D. pulex, one third of its genome comprises lineage-specific genes that lack orthologs in other eukaryote genomes and thus lack any functional annotation [20], and the situation is probably similar for the D. magna genome [17]. Transcriptomic studies in D. magna and related species have been up to now mostly conducted to evaluate the response of genes to environmental factors with few attempts to address transcriptomic changes across life-cycle stages [1,[25][26][27][28][29][30][31][32][33][34]. This means that little is known about the role of genes during its development, reproduction, and across sexes, which may compromise ecological and ecotoxicological interpretations.
D. magna has a well-known development and lifecycle stages, which combine parthenogenetic and sexual reproduction [35]. Parthenogenetic reproductive females allocate freshly formed eggs into the brood pouch every 3-4 days at 20°C, under optimal food conditions [10]. Freshly released eggs develop independently of the mother across 12 visually recognizable embryonic stages into free swimming juveniles [36] that are released by the mother just before it moults and allocates a new clutch of eggs into its brood pouch [10]. Juveniles, which are morphologically similar to adults, moult four or more times, depending on food conditions, before closing the circle by allocating the first clutch of eggs into their brood pouch [10]. In the ovary of a D. magna reproductive female co-exist clusters of pre and vitellogenic oocytes [35], the latter being provisioned with yolk during the intermoult interval to form the new clutch of eggs at the end of it, whereas the former will mature and became vitellogenic in the next reproductive cycle. Upon adverse conditions, D. magna females produce males whose haploid sperm can fecundate haploid oocytes to form an ephippia containing two sexually produced eggs, a process that can be triggered independently from environmental conditions by exposure to the crustacean juvenile hormone methyl farnesoate [35,37]. Ephippia eggs can survive for decades in water bodies' sediments and develop into free-swimming female juveniles when conditions are favourable. Thus D. magna is also an excellent model to study large-scale gene transcription across development, reproduction and sexes to provide global relationships between transcription and life-cycle stages.
Novel developments of bioinformatic tools [38] and of DNA/RNA-seq technologies have allowed the publication of the Daphnia pulex full genome [20] and more recently of the full D. magna transcriptome [1]. Information reported by Orsini et al. [1] of approximately 42,000 gene validated models, including splicing forms, was used in this present study to design a new microarray, representing the full gene set of the D. magna transcriptome coupled with most up-to-date functional information. The new microarray platform was then used to study gene transcription patterns for the full transcriptome of D. magna across early and late embryonic stages, juveniles, reproductive females and males (Fig. 1).

Microarray performance
The probe design was checked for specificity performing a complete probe match analysis against the D. magna original gene models. Using exact match criteria 3616 probes match to two or more transcripts. If allowing for 1-3 mismatches (5% divergence) 5371 of the probes matched to two or more transcripts.
The microarray ability of representing real changes in transcription was checked performing qPCR on four relevant genes from three different pathways (Krebs cycle, tryptophan metabolism and ecdysteroid pathways). The results between the two technological platforms were highly correlated (Pearson correlation > 0.73, P < 0. 01, N = 20, Additional file 1: Figure S1).

Differentially transcribed genes
About 109,640 probes showed fluorescence values above the background. These probes belonged to 38,449 unique genes out of the 41,317 D. magna genes present in the microarray design (see Additional file 2 for further information regarding the non-expressed genes). Hierarchical clustering of fluorescence values for the 109,640 probes in all samples showed a good grouping of biological replicates (Fig. 2), a sign for good quality of the design. The analysis defined three main clusters: Cluster 1 composed by F2, F3 samples; Cluster 2 grouping E3, M, F 1 and J stages; and Cluster 3 including only E 1 samples.
The total number of differentially transcribed genes (DEG) can also be an important measure of how active the cells of an organism are during a particular life stage. The number of probes (DEP), and the unique genes these represent, differentially transcribed across lifestages was quite variable ( Fig. 3 and Additional file 1: Tables S1 and S2). The highest number of DEP corresponded to the E1 stage (17881 DEP, representing 9663 unique genes), steadily decreasing to 10467 DEP (6013 Genes) in E3 and 4332 DEP (2893 Gene) in juvenile stages (Fig. 3). F1 females had the lowest number of DE compared to the all other life stages with 2100 DEP (1294 Genes), being that number increased to 3439 DEP (2292 genes) and 4141 (2900 genes) in F2 and F3, respectively. To allow a better comparison of transcriptomic changes across sexes all female stage x replicates were considered together, which totalled 14272 DEP (8026 genes). Furthermore, up to 6352 genes (11915 DEP) were differentially upregulated in both E1 and in F3 (E1F3, Fig. 3). Finally, males differentially expressed a total of 3980 DEP (2983 genes). In most life stages there were more up-regulated DEG than down regulated ones except for F1. Details of the expression levels of differentially transcribed probes and genes at the various studied life-stages are provided in Additional file 3.
The heat map of DEG included three clusters (Fig. 4a): the first one constituted by E1 samples, a second other including E3 and M samples, and a third cluster grouping the three female stages (J, F1, F2 and F3). Transcription patterns of up and down regulated DEGs in E1, E3, J, F1, F2, F3, F123 and M were stage-specific, whereas that of E1F3UP were up regulated in E1 and F3 (E1F3) and to a lesser extent in F2 and F1.

Functional analysis
A total of 25,019 (60%) genes present in the microarray were functionally annotated to some level by using the three major databases (Drosophila, REFSEQ and Swissprot KB), being 16,081 genes mapped to at least one GO Term and 7995 annotated to a KEGG-curated gene. Using a high stringency probability (p < 0.001), a total of 2905 distinct GO terms were enriched across the b a Fig. 1 Diagram representing the experimental design. a Seven life-stages were selected: embryos < 12 h after oviposition (Eggs 1, E1; corresponding to eggs of stage 1 following [1]), embryos within the last 12 h to be release as free swimming first instar juveniles (Eggs 3, E3; corresponding to embryos of stage 12 according to [1] (Table 1). Selected GO terms and their coverage in each stage are shown in Table 2.
KEGG mapping of specific DEGs generally agreed with the identified GO terms, but they present a better overview on how the number of differentially expressed genes mapped to specific pathways along the development of the organism ( Fig. 5 and Additional file 5).

Discussion
Most of the probes included in the microarray match a unique gene transcript with only about 3.2% of probes matching two or more gene transcripts. Furthermore correlation values of selected gene transcript levels between qPCR and microarray probes were high, even for low differentially transcribed genes, thus indicating a good quality of the probes, capable of reflecting even small changes in mRNA abundance.
More than 93% of the genes that are reported to be transcribed in D. magna genome [1] showed fluorescence values above background levels in at least one out of the seven studied life-cycle stages. Hierarchical clustering of samples and the number of differentially transcribed genes (Figs. 2 and 3) separated eggs from stage E1 that showed the greatest levels of differential transcribed genes (34.5% of the 38,449 genes), reproductive females of stage F2 and F3 with up to 18.5% of the transcribed genes and the rest of life stages (E3, J, F1 and M) with the remaining 47% of the gene transcripts. These three clusters had distinctive developmental/physiological characteristics. D. magna females transfer large amounts of resources including mRNA to their eggs to ensure their development and survival upon hatching [28,[39][40][41][42]. Therefore, E1 is likely to be enriched with genes associated to embryonic development and maternal transfer. Indeed up to 6352 genes were upregulated in E1 and in females F3 (E1F3 group in Fig. 3). This is in line with reported transcriptomic results in D.
melanogaster, Anopheles and C. elegans, which found that the majority of genes were transcribed during early embryogenesis stages due to the contribution of both maternal and zygote ones [2][3][4][5]. The second cluster, grouping F2 and F3 females, can be considered representative of the most reproductive active stages, whereas life-stages included the third cluster (E3, J, F1 and M) are associated to high moulting activity and low-reproduction related processes [43,44]. The relative low levels of differentially transcribed genes in F1 females (4.6% relative to the total), combined with an overrepresentation of down-regulated genes (see Additional file 1: Table S2) can be related to the fact that F1 organisms perform a reset of their metabolic and physiological processes at the beginning of the intermoult cycle [43,44]. Therefore, most of the specific differences showed in this stage could be related to de-activation of gene transcription.
Gene transcription across the studied life-stages ( Fig. 4) also evidenced clear distinctive patterns across egg, juvenile and adult stages. Up regulated DEGs in reproductive females and males showed low transcription levels in the rest of life-stages whereas low transcribed ones were highly transcribed in the remaining stages. In the embryonic stage E3 and juveniles, genes having high or low levels of transcription changed little in the remaining stages. The early embryonic stage E1 showed two distinctive transcription patterns: a cluster of 9663 unique genes highly or poorly transcribed in that stage and a cluster of 6352 genes highly transcribed in that stage and in female stages F3 and F2. These results agree with those reported in Drosophila and C. elegans supporting the argument that many transcripts in E1 have a maternal origin [2,3]. The unique gene patterns identified in most life-stages and across sexes also support previous studies that associate these stages with specific subsets of genes being expressed in specific tissues and biological processes [2,3].

Gene function analysis
E1 transcriptome included the highest number of significant GO terms (1380), mainly in categories related to division and proliferation, morphogenesis, homeostasis, and immune response ( Table 2). Homeostasis is considered important to cells on high metabolic rates in which accumulation of some metabolites can become toxic. The enrichment in methylation-related transcripts may be linked to epigenetic processes, which are common during embryo development in Daphnia and other arthropods [2,28,45]. Interestingly the cluster of genes from E1 coming from a maternal origin (E1F3) were enriched with 829 GO terms, 50% of them also present in E1 (Tables 1 and 2). Transcripts enriched in E3, with up to 1323 significant GO terms, also had important representation of GOs related with morphogenesis/organogenesis. They also include several GOs related to  response to chemicals (Table 2), perhaps preparing the organism to survive in the outside environment [46]. This stage also show an enrichment of GO terms related with sensorial and locomotion processes (neuronal activation and differentiation, hormonal metabolic processes and binding and locomotor behaviour), which would prepare the embryo for the free-swimming stage. The E3 transcriptome was also enriched in GOs related to cuticle constituents and elasticity, a functional category common to all the adult female stages. Transcriptome from the juvenile stage was annotated to 97 significantly enriched GO terms (Additional file 4), including the distinctive GO term "Bursicon neuropeptide hormone complex" (Table 2), which is related to the hardening of carapace [47]. This is in line with the fact that moulting and growth processes dominate in this stage over other biological functions such as reproduction [48]. Enriched GO terms from F1 transcriptome (145 GO terms in total) included two distinctive GO's: dormancy process and activin complex. The first GO term is possibly related to the re-setting of moulting and reproductive processes in which most bodily functions are "on hold" [49]. The activin complex GO term is directly related with the Follicle Stimulating Hormone (FSH) a key hormone to trigger the development of ovules [50]. This is a key process of extreme importance to the reproduction of these organisms showing that any impairment of these genes during this life stage period can be critical for the reproductive output. This is consistent with the fact that ovulation occurs precisely at the F1 stage [35]. In F2, specific transcripts were annotated to 146 significant GO terms, with one distinctive term, hair cycle phase, which can be linked in these organism to the formation of carapace and, possibly, also to hardening processes. Analysis of the F3 transcriptome identified 91 enriched GO terms, two of them with a high coverage (47%, 71% coverage of the annotated probes DEG in this stage) and belonging to "Anatomical Structure Development or morphogenesis". These terms may be related to two processes: (1) the development of the ovary and formation and provisioning of eggs and (2) the maternal transfer of RNAs from the mother to the newly formed eggs to start development. Interestingly, a large fraction of probes annotated to Egg stages and of maternal origin (E1F3) also belong to this GO term. Another key term uniquely related to F3 was "Abcission", which in arthropods is involved in carapace shedding [51], a key event in the development of crustaceans that are only able to grow by periodically releasing the old carapace and the subsequent formation of a new one. This last GO term is likely related to the fact that females F3 are close to a moult event. Differentially transcribed genes of reproductive females (F123) belonged to 1004 GOs. Almost half (444) of those GO also existed in E1 and E3. Malespecific transcripts were annotated to 97 GO terms, including the unique category of "Male Germ Cell proliferation", which in turn includes the gene doublesex, linked in Daphnia to male sexual differentiation [52]. Males also present a high number of probes associated with carbohydrate transport activity, like F3, which can related to a high demand of energy for swimming in males [53] and to reproductive activities for females [48]. The metabolism-oriented KEGG pathway mapping ( Fig. 5 and Additional file 5) shows a high number of DEG related to carbon, starch and sucrose, tryptophan and fatty acid metabolism in embryo (E1, E3 and DEGs related to DNA replication were abundant in all life stages and in the cluster of E1 from maternal origin (E1F3) with except for F2. In the parthenogenicallyreproducing sister species Daphnia pulex, RNA and DNA concentrations were related to changes in metabolic activity associated with moulting cycle and ontogenetic development [54]. Ontogenic developmental processes are likely to dominate in embryo stages (E1, E3). Females F2 are just in the middle of the intermoult cycle, whereas the rest of free living stages (J, F1, F3) are likely to be just at the beginning or end of the moulting cycle.
An important point to consider when evaluating the quality of the data we are reporting of its functional analysis is the fact that males have no DEG mapped on females processes and hormones (Gonadotropin, progesterone, oocyte and ovarian steroidogenesis). All these processes were highly enriched in eggs, juveniles and F1,2,3.
Embryonic stages showed a high expression of genes annotated to many signalling pathways, like MAPK (response to stress), RAS (cellular signal transduction), RAP (vital for effective signal transduction), cAMP (second messenger important in several biological processes), and Wnt (signal transduction pathways made of proteins, extremely important for embryonic development). Genetic and molecular studies in Drosophila melanogaster over the last several years show that these signalling pathways are functionally conserved and that they participate in numerous processes during normal development [55][56][57][58][59].
Similarly, genes implicated in neuronal processes were highly expressed in egg and female stages. This is in line with the development of the neuronal system in metazoans along the embryonic stages [60]. Finally, hormonal processes related to Prolactin and insulin were evenly differentiated across life cycle stages, whereas genes belonging to thyroid hormone synthesis and oxytocin KEGG pathway, were differentially regulated in egg and female stages. Despite that genes and peptides of the above signalling pathways have been found in D. pulex and D. magna genome, little is known about their expression and function [61][62][63].

Conclusions
In this work we present a new microarray platform designed using the full set of gene models representing the complete genome of D. magna [1]. Up to 93% of the existing 41,317 D. magna gene models showed differential transcription patterns across the developmental stages of D. magna, from which more than 59% were functionally annotated. The transcriptome analysis show a large number of embryo-specific transcripts, likely related to the high metabolic rates needed for cellular multiplication and morphogenesis of several body organs, as inferred by the high number of DEG related to DNA, RNA and ribosome biogenesis. In adult females, most DEG were involved in reproductive processes and carapace shedding. Differentially transcribed genes in males were enriched in specific genes of male sexual differentiation genes, like doublesex [52]; consistently, they did not include hormonal female-related functional categories. These results are in line with those reported in Drosophila, showing that about 90% of transcription of embryo-specific transcripts occurs during the first hours of embryo development [2].

Experimental animals and culture conditions
All experiments were performed using a wellcharacterized single clone of D. magna (Clone F), maintained indefinitely as pure parthenogenetic cultures. Individual or bulk cultures of 10 animals/L were maintained in ASTM hard synthetic water [64] as described in [44]. Individual or bulk cultures were fed daily with Chorella vulgaris Beijerinck (5 × 10 5 cells/mL, respectively, corresponding to 1.8 μg C/mL; [10]). C.
vulgaris was grown axenically in Jaworski/Euglena gracilis 1: 1 medium (CCAP, 1989). Algae were harvested in the exponential phase of growth, centrifuged and then re-suspended in ASTM hard water. The number of algal cells in freshly prepared medium was checked daily from absorbance measurements at λ = 650 nm in a dual-beam spectrophotometer (Uvikon 941) using standard calibration curves based on at least 20 data points, with an r 2 > 0.98. The culture medium was changed every day, and neonates were removed within 24 h. Photoperiod was set to 14 h light: 10 h dark cycle and temperature at 20 ± 1°C.
For this study, females across six different developmental stages and one adult male stage were collected in triplicate. Embryonic samples were taken according to the embryonic development staging works by Mittmann et al. [36], while the free swimming stages were sampled following major physiological events. A graphic representation of these stages is depicted in Fig. 1. The first egg stage (E1) includes embryos < 6-8 h after ovideposition in stage 1 [36]. At Stage 1 the nucleus lies within non-transparent globular pale to greenish yolk granules and several larger oil drops. The egg is covered tightly by the chorion. The second embryonic stage selected (E3) corresponds roughly to stage 12 [36] and includes embryos within their last 12 h before the release from the mothers brood pouch. At this stage the embryonic development is mostly completed and is in the path of transformation to the first instar juvenile [36], being already able to swim. The choice of this stage was made to represent a neonate without the interference of feeding behaviour and or digestion processes. The third selected stage (J) is a fourth instar juvenile or pre-adolescent stage. Animals of about 4 days old (with no ovary visible) [65], targeted to represent a full working organism but without the interference of reproductive processes.
Adult third clutch female stages (F1), (F2) and (F3) represent, respectively, reproductive females at the beginning < 12 h), middle (30 h < X < 40 h) and end (> 66 h) of their reproductive/intermoult cycle, which usually lasts 3 days at 20°C with non-limiting food conditions as in the present study [44]. Females F1 were collected just after moult and release of the new clutch of newly formed eggs of stage 1 (i.e. E1 embryos) into their brood pouch (< 8 h). Females F2 were sampled after approximately 36 h and F3 were collected within the last 24 h before moulting and releasing fully developed embryos. Females F3 thus hold embryos of stage 12 (our E3 embryos). Ovulation and egg provisioning occurs simultaneously during a single reproductive/intermoult cycle [35]. In females F1 pre-vitellogenic eggs pass to vitellogenic ones in the ovary and starts the egg provisioning, which is fully accomplished in F3. The development of embryos in the brood pouch also occurs simultaneously as ovulation and egg provisioning [35].
Females were de-brooded before being collected, by gently flushing water into the brood pouch. Removing the developing embryos from the brood pouch also removes their RNA, allowing for a much more robust profiling of the female adult stages [10]. The final samples include reproductively active males. These males were obtained exposing parthenogenic females to 0.5 μM of methyl farnesoate [37]. After birth, the males were selected observing the size of the first antenna and let to grow until 14 days, representing a physiological status roughly coinciding with female stages.

RNA extraction
Total RNA was isolated from the samples using Trizol reagent (Invitrogen, USA) and following manufacturer protocols with slight modifications. After RNA isolation, DNAse treatment was performed according to manufacturer protocols, followed by a double phenol-chloroform and another chloroform extraction for further purification. RNA was precipitated using Sodium acetate and 100% ethanol, re-suspended in RNAse free water being quantity and quality checked using a NanoDrop D-1000 Spectrophotometer (NanoDrop Technologies, USA). Samples presenting a ratio 230/260-260/280 between 1. 9-2.1 were selected. RNA integrity was checked using an Agilent 2100 Bioanalyzer (Agilent Technologies, USA). Only the samples showing invertebrate-adjusted RIN values above 9 were used for microarray analysis.

Microarrays
This new design was made using the gene models developed by Orsini et al. [1]. This gene set includes 41,317 gene models representing the full transcriptome of Daphnia magna. Four probes were designed for every gene, using the algorithms available in the e-array tool from Agilent (Palo Alto, USA). One of the probes was the best possible probe, while the other three were designed to represent a best geographical representation of the gene sequence, i. e., one probe representing the beginning of the sequence, another one representing the central part and the last the end part of the sequence. After probe design, redundancy was checked and in case of existence, redundant ones were removed. A set of 3500 random computer-generated probes were also added as further negative controls. Further e-array based quality controls were added, resulting in a microarray with 165,000 probes. This was then printed on a 4 × 180,000 format (Agilent 66414 design; GPL22721). A total of three replicates per treatment were used. One μg of total RNA was used for all hybridizations. cDNA synthesis, cRNA labelling, amplification, and hybridizations were performed following the manufacturer's kits and protocols (Quick Amp labelling kit; Agilent, Palo Alto, CA). The Agilent one-color Microarray Based Gene Expression Analysis v6.5 was used for microarray hybridizations according to the manufacturer's recommendations. Microarray images were generated by an Agilent high-resolution C microarray scanner. Data was resolved from microarray images using Agilent Feature Extraction software v10.7. Raw microarray data from this study have been deposited at the Gene Expression Omnibus Web site (www.ncbi.nlm. nih.gov/geo/) with accession number GSE90810.

Gene expression analysis
Microarray data were analysed using Gene Spring GX v13.0 software (Agilent, USA). Fluorescence data were normalized using quantile normalization and baseline transformation to the median of all samples.
The quantile 95 of the 3500 negative probes was calculated and this value was assumed as being the fluorescence background noise value of each sample. Only the probes having fluorescence values above background noise values (109,640 probes in total) were considered for further analyses. Sample Clustering of the total expressed genes across life-stages was analysed using the Multi-Experiment viewer MeV4 software [66] by hierarchical clustering using Pearson correlation algorithm.
Differentially transcribed probes (DEP) of D. magna individuals in a given life-stage were identified by pairwise comparison of the normalized replicated fluorescent levels of that stage against the other stages using Student t test p < 0.001) and a fluorescence change cut off of 1.5-fold. The later term was calculated as the quotient between the power of two of the mean normalized fluorescence of both stages involved in the comparison. Only those probes whose normalized fluorescence changed significantly ≥1.5 fold in all pairwise "stage" vs "all others" comparisons were considered. For a given stage probes having greater or lower normalized fluorescence values versus the rest of stages were considered up and down regulated. This approach was found to be more conservative than the combination of ANOVA approaches and the Benjamini-Hochberg false discovery rates correction. Finally, differentially transcribed genes (DEG) was defined as a non-redundant representation of the genes uniquely represented by at least one DEP. To determine genes present in E1 from maternal origin, differentially transcribes probes upregulated in both E1 and F3 were determined by comparing the normalized replicated fluorescent levels of both stages against the juvenile stage using Student t test p < 0.001. Stages F3 and J were selected since F3 is the female stage where more resources from ovaries are allocated to eggs, whereas J can be considered an immature female stage [10].

Validation of microarray results by qPCR
Microarray results were validated with real-time quantitative polymerase chain reaction (qPCR). We selected four differentially expressed genes from different pathways/gene families which had already been successfully used in previous studies [46,67]: Krebs cycle (isocytrate dehydrogenase, idh; ATP citrate lyase, ATPCL), tryptophan metabolism (dopamine decarboxylase, Ddc) and Ecdysteroid pathways (Ecdysone receptor, EcR b). The G3PDH gene (glyceraldehyde 3-phosphate dehydrogenase) was used as internal control (Housekeeping). Primers for each one of these genes were designed with Primer Express® Software v3.0.1(Thermofisher, USA) and are provided elsewhere [46,67]. qPCR was performed according to manufacturer's protocols.

Gene model annotation
For the optimal interpretation of transcriptional changes it is fundamental to have a functional annotation of all genes. To achieve this, all the 41,317 genes were searched for homology and functional annotation, including Gene Ontology (GO) and Kyoto Encyclopaedia of Genes and Genomes (KEGG). This was performed using Blast2Go Suit [26]. To maximize homology findings, this search was performed using BLASTX alignment algorithms (search protein databases against a translated nucleotide query). Search was performed using as target the three main databases of interest to our platform: Drosophila, RefSEQ and SwissprotKB. The three annotation sets were finally merged. The ultimate objective of re-annotation of the gene models was to improve as much as possible the functional annotation of these. This is the reason why we limited our search to the reference databases as these are the ones with most functional mapping. The full details of the annotation are available in GEO together with the microarray design (GPL22721).

Additional files
Additional file 1: Word file with the qPCR validation of the results in Figure S1 and the total number probes and genes of a given life-stage differential transcribed across the remaining ones and those up and down regulated in Tables S1 and S2 (DOCX 40 kb)