Taenia

A custom microarray was designed from the sequence of the junco transcriptome. For 33,545 contigs, three unique probes are present on the array, while another 61 contigs are represented by two probes and 65 contigs are represented by one probe, accounting for 100,822 probes on the array. An additional 34,365 probes were selected from the remaining singletons (one probe per chosen singleton). The array also contains control probes and 2,604 random probes designed to reflect the genome nucleotide composition by Markov modeling to experimentally determine the appropriate thresholds that measure significant hybridization signals over the background. Thus, each sub-array consists of over 137,000 long-oligonucleotide (60 bp) probes, and 12 such sub-arrays are placed on each glass slide (Roche NimbleGen Inc., Madison, WI). The microarray platform is deposited at NCBI Gene Expression Omnibus (GEO; accession number GPL14995).
We collected adult dark-eyed juncos from breeding grounds near Mountain Lake Biological Station (Pembroke, VA) in mist-nets between May 7 and 14, 2010 and held them individually in a semi-naturalistic outdoor aviary where they were neither acoustically nor visually isolated from other juncos, as part of a larger experiment. On June 9 and 10 individuals were euthanized by overdose of isoflurane. Tissues, including whole brains, were collected rapidly and stored on powdered dry ice within 20 min post-mortem to ensure negligible RNA degradation [62 ]. Brains were dissected into 14 distinct regions using anatomical landmarks, following previously established methods [63 (link)] based on the zebra finch brain atlas. These brain regions included the hypothalamus and the ventral medial telencephalon (VmT), which primarily consists of the nucleus taeniae, the avian homologue of the medial amygdala [64 -66 (link)].
RNA from VmT, hypothalamus, liver, and pectoralis was extracted in TRIzol® following manufacturer's directions (Invitrogen, Carlsbad, CA). The microarray protocol follows previously published methods [67 ]. Briefly, total RNA was reverse-transcribed to ss-cDNA in the presence of oligodT primer and SuperScript II reverse transcriptase. This ss-cDNA was then converted to ds-RNA and labeled using CY-labeled random nonmer primer (either Cy3 or Cy5) and Klenow fragment (following NimbleGen labeling protocols). We then hybridized 4 g of each of two labeled samples (one Cy3, one Cy5) to each sub-array and followed manufacturer's directions for post-hybridization washing and scanning (Roche NimbleGen, Inc., Madison, WI). Imaging was accomplished by Axon GenePix 4200A scanner (Molecular Devices, Sunnyvale CA) with GenePix 6.0 software and data were extracted with NimbleScan 2.4 (Roche NimbleGen, Inc., Madison WI). Raw microarray data were processed with the limma package [68 ] in R version 2.9.0 [69 ] to normalize expression scores.
To determine if a gene was expressed, we calculated the 97.5% quantile for expression score of random probes in each individual as the cutoff for calling expression. Thus, for each individual, a called expression is significant at a p-value of 0.025. For each contig, we tested the median probe value against this threshold, and for singletons we used the single expression value. Because our design employed biological replicates, we called a contig or singleton expressed only if at least three of the six individuals in a group were called as expressed, thus reducing the p-value further to 0.0006 (the probability of obtaining at least three of six individuals called for expression of a random probe). From this, we determined whether or not a gene had expression support in any of our tissues-sex pairings, and which genes were restricted to expression in one sex.

Two phylogenetic analyses were performed using putative OGPs obtained from the genomes (g), transcriptomes (t), or EST sequences of 11 species of Neodermata: the monopisthocotyleans Rhabdosynochus viridisi (t), Scutogyrus longicornis (t), Gyrodactylus salaris (g), and Neobenedenia melleni (EST) [28 (link),29 (link)]; the polyopisthocotyleans Protopolystoma xenopodis (g) and Eudiplozoon nipponicum (t) [29 (link),30 (link)]; the cestodes Echinococcus multilocularis (g), Hymenolepis microstoma (g), and Taenia asiatica (g) [29 (link)]; and the trematodes F. hepatica (g) and Schistosoma mansoni (g) [29 (link)]. The free-living platyhelminths Bothrioplana semperi (g), belonging to the order Bothrioplanida, and Schmidtea mediterranea (g), belonging to the order Tricladida, were used as outgroups [12 (link),29 (link)]. The first analysis included single-copy OGPs retrieved from BUSCO v4 [31 (link)], using the core metazoan dataset, which contains 978 genes and the script BUSCO_phylogenomics.py with parameters “-supermatrix” and “-psc 70”. The second analysis included simply OGPs (OMA Groups) obtained through OMA Standalone [32 (link)], using the script filter_groups.py [19 (link)]. OMA Groups contained a maximum of one representative gene per species. When multiple co-orthologs exist, OMA selected one to be in the OMA Group. Only the OGPs present in at least 10 species were included in the analyses.
The BUSCO proteins used in the phylogenetic analyses were annotated with the odb10 database [33 (link)], whereas the OMA proteins were annotated using BLASTp [34 (link)] against the UniProtKB/Swiss-Prot database (e-value < 1 × 10⁻⁴). In addition, the proteins were mapped to Gene Ontology (GO) terms using the PANNZER2 web server [35 (link)]. The visualisation and GO term enrichment analysis were performed in WEGO [36 (link)].
The OGPs were aligned with Muscle v3.8.31 [37 (link)], trimmed with trimAL [38 (link)] using the automated mode (-automated1), and concatenated. The best evolutionary models were obtained with the ModelFinder program [39 (link)]. Phylogenetic trees were constructed in IQ-TREE v1.6.12 [40 (link)], using the Shimodaira–Hasegawa-like approximate likelihood ratio test (SH-aLRT) (1000 replicates) and 1000 ultrafast bootstrap approximations to calculate the support values of the clades. Phylogenetic trees were also constructed in RAxML v8 [41 (link)], with 1000 bootstrap (Bs) iterations to calculate the support values of the clades. The trees were visualised with FigTree v1.4.2 (http://tree.bio.ed.ac.uk/software/figtree/, accessed on 1 November 2021). The trees were constructed using the models LG + F + R4 and LG + F + G + I for IQ-TREE and RAxML, respectively. Constrained trees were constructed in IQ-TREE (using the -g option), using alternative scenarios obtained from the literature as a guide (see Table 1) and the LG + F + R4 model. To determine whether the topology found in this study was significantly better than constrained trees, tree topology tests were performed in IQ-TREE using the LG + F + R4 evolutionary model with the options: -au, -zb.

18S rDNA sequences of Dipylidium caninum, Hymenolepis nana, Taenia solium, Taenia saginata, Taenia taeniaeformis, Echinococcus granulosus, Echinococcus multilocularis, Fasciola hepatica, Fascioloides magna, Ascaris lumbricoides, Ascaris suum, Enterobius vermicularis, Toxocara canis, Toxocara cati, Trichuris trichiura, Trichuris vulpis were used to design a six-primer set. Four 18S rDNA sequences of Giardia intestinalis were used to design PCR primers. Detailed data on the parasites and 18S rDNA sequences are presented in S4 Table. Partial regions (250bp) of 18S rDNA of all tested helminths were amplified using the following sets of primers: CF1 (5’ GCGGGGRCGTTTGTATGGCTGC 3’), NF1 (5’ ACGGGGRCATTCGTATyGCTGC 3’), TtF 1 (5’ GCGGGGACGTTTGTATGGTTGCG 3’), FhF 1 (5’ ACGGGGGCATTTGTATGGCGGT 3’), CNR1 (5’ CAACCATACTTCCCCCGGAACCSAAA 3’) and TR1 (5’ CCATACTTCCCCCGGAGCCCAAA 3’). PCR was performed in a TM 100 (BioRad) thermal cycler. The reactions were conducted in a 50 μl reaction mixture containing 5.0 μl of DNA template, 1 μl (1U) of Color Taq DNA Polymerase (EURx), 1 μl of dNTPs mix (10 mM), 0.5 μl of TR1, TtF1 and FhF1, 1 μl of CF1 and CNR1, 1.5 μl of NF1 primer (20 mM), 5 μl of 10 × Polymerase buffer (pH 8.6, 25 mM MgCl ₂) and 33.0 μl of MiliQ water. A negative control—nuclease-free water was added to the PCR mix instead of the tested DNA. DNA amplification was performed according to the following program: denaturation at 95°C for 1 min, followed by 34 cycles of denaturation at 95°C for 10 s, annealing at 58°C for 10 s and extension at 72°C for 20 s, with a final extension performed at 72°C for 3 min.