Rhea

The Rhea pipeline consists of six main R scripts that perform the tasks of normalization of input tables, calculation of alpha-diversity, beta-diversity, taxonomic relative abundances, serial group comparisons, and correlations. These scripts rely on the R packages ade4, GUniFrac, phangorn, Hmisc, corrplot, plotrix, PerformanceAnalytics, reshape, ggplot2, gridExtra, grid, ggrepel, gtable, and Matrix that by themselves have several dependencies. The installation of the packages is performed automatically within the scripts when run for the first time. For the purpose of demonstrating and illustrating the different features of Rhea, the publicly available sequence data from the study by Müller et al. (2016) (link), with ENA accession PRJEB13041, were analysed with the web platform IMNGS (http://www.imngs.org) and the output OTU-table and files (also available for download through the GitHub repository) were used for analysis. In this template study, the impact of housing conditions and diet on mouse faecal microbiota and gut barrier was investigated. For demonstration of the variability of OTU-specific relative abundances among technical replicates, a small study of 10 amplicon libraries constructed from the same human faecal sample were sequenced by Illumina MiSeq. The raw data of this study were deposited to ENA and are available under accession PRJEB14963. Those data were also processed with the web platform IMNGS and the OTU-table and all intermediate steps for the analysis of Coefficient of Variation can be find in Table S1.
The following sections describe in detail each of the main functions of Rhea, thereby emphasizing on important concepts underlying data processing using the scripts. Meticulous documentation of all scripts is provided online at the link given above in the abstract. To minimize manual handling of data, intermediate files generated during processing are automatically transferred to folders where they are needed for downstream analysis, on the condition that the original folder structure of Rhea is kept unchanged. Illustrations shown in the present manuscript correspond to raw outputs generated by Rhea, with very minor post-production manual changes (i.e., size and orientation may have been changed to facilitate publishing).

To test for universality of primers and cycling conditions, we performed parallel experiments in three different laboratories (Berkeley, Cologne, Konstanz) using the same primers but different biochemical products and thermocyclers, and slightly different protocols.
The selected primers for 16S [30 ] amplify a fragment of ca. 550 bp (in amphibians) that has been used in many phylogenetic and phylogeographic studies in this and other vertebrate classes: 16SA-L, 5' - CGC CTG TTT ATC AAA AAC AT - 3'; 16SB-H, 5' - CCG GTC TGA ACT CAG ATC ACG T - 3'.
For COI we tested (1) three primers designed for birds [7 (link)] that amplify a 749 bp region near the 5'-terminus of this gene: BirdF1, 5' - TTC TCC AAC CAC AAA GAC ATT GGC AC - 3', BirdR1, 5' - ACG TGG GAG ATA ATT CCA AAT CCT G - 3', and BirdR2, 5' - ACT ACA TGT GAG ATG ATT CCG AAT CCA G - 3'; and (2) one pair of primers designed for arthropods [2 (link)] that amplify a 658 bp fragment in the same region: LCO1490, 5' - GGT CAA CAA ATC ATA AAG ATA TTG G - 3', and HCO2198, 5'-TAA ACT TCA GGG TGA CCA AAA AAT CA-3'. Sequences of additional primers for COI that had performed well in mammals and fishes were kindly made available by P. D. N. Hebert (personal communication in 2004) and these primers yielded similar results (not shown).
The optimal annealing temperatures for the COI primers were determined using a gradient thermocycler and were found to be 49–50°C; the 16S annealing temperature was 55°C. Successfully amplified fragments were sequenced using various automated sequencers and deposited in Genbank. Accession numbers for the complete data set of adult mantellid sequences used for the assessment of intra- and interspecific divergences (e.g. in Fig. 5) are AY847959–AY848683. Accession numbers of the obtained COI sequences are AY883978–AY883995.
Nucleotide variability was scored using the software DNAsp [31 (link)] at COI and 16S priming sites of the following complete mitochondrial genomes of nine amphibians and 59 other vertebrates: Cephalochordata: AF098298, Branchiostoma. Myxiniformes: AJ404477, Myxine. Petromyzontiformes: U11880, Petromyzon. Chondrichthyes: AJ310140, Chimaera; AF106038, Raja; Y16067, Scyliorhinus; Y18134, Squalus. Actinopterygii: AY442347, Amia; AB038556, Anguilla; AB034824, Coregonus; M91245, Crossostoma; AP002944, Gasterosteus; AB047553, Plecoglossus; U62532, Polypterus; U12143, Salmo. Dipnoi: L42813, Protopterus. Coelacanthiformes: U82228, Latimeria. Amphibia, Gymnophiona: AF154051, Typhlonectes. Amphibia, Urodela: AJ584639, Ambystoma; AJ492192, Andrias; AF154053, Mertensiella; AJ419960, Ranodon. Amphibia, Anura: AB127977, Buergeria; NC_005794, Bufo; AY158705; Fejervarya; AB043889, Rana; M10217, Xenopus. Testudines: AF069423, NC_000886, Chelonia; Chrysemys; AF366350, Dogania; AY687385, Pelodiscus; AF039066, Pelomedusa. Squamata: NC_005958, Abronia; AB079613, Cordylus; AB008539, Dinodon; AJ278511, Iguana; AB079597, Leptotyphlops; AB079242, Sceloporus; AB080274, Shinisaurus. Crocodilia: AJ404872, Caiman. Aves: AF363031, Anser; AY074885, Arenaria; AF090337, Aythya; AF380305, Buteo; AB026818, Ciconia; AF362763, Eudyptula; AF090338, Falco; AY235571, Gallus; AY074886, Haematopus; AF090339, Rhea; Y12025, Struthio. Mammalia: X83427, Ornithorhynchus; Y10524, Macropus; AJ304826, Vombatus; AF061340, Artibeus; U96639, Canis; AJ222767, Cavia ; AY075116, Dugong; AB099484, Echinops; Y19184, Lama; AJ224821, Loxodonta; AB042432, Mus; AJ001562, Myoxus; AJ001588, Oryctolagus; AF321050, Pteropus; AB061527, Sorex; AF348159, Tarsius; AF217811, Tupaia; AF303111, Ursus (for species names, see Genbank under the respective accession numbers).
16S sequences of a large sample of Madagascan frogs were used to build a database in Bioedit [32 ]. Tadpole sequences were compared with this database using local BLAST searches [33 (link)] as implemented in Bioedit.
The performance of COI and 16S in assigning taxa to inclusive major clades was tested based on gene fragments homologous to those amplified by the primers used herein (see above), extracted from the complete mitochondrial sequences of 68 vertebrate taxa. Sequences were aligned in Sequence Navigator (Applied Biosystems) by a Clustal algorithm with a gap penalty of 50, a gap extend penalty of 10 and a setting of the ktup parameter at 2. PAUP* [34 ] was used with the neighbor-joining algorithm and LogDet distances and excluding pairwise comparisons for gapped sites. We chose these simple phenetic methods instead of maximum likelihood or maximum parsimony approaches because they are computationally more demanding and because the aim of DNA barcoding is a robust and fast identification of taxa rather than an accurate determination of their phylogenetic relationships.

Details on the generation of new rhea (talin) and Vinculin alleles can be found in Supplemental Experimental Procedures.
For wing blister quantification, mitotic clones were generated in the wings of heterozygous flies by crossing rhea mutant males to w; P{w[+], Gal4}Vg[BE] P{w[+], UAS::FLP}; P{FRT}2A (with the white+ excised from P{FRT2Aw[hs]}) females. Embryonic phenotype quantification was performed on mutant embryos lacking both maternal and zygotic wild-type talin and/or vinculin, as they were obtained from germline clones generated in heterozygous mutant females by crossing rhea mutant females (with wild-type Vinculin or ΔVinc) to P{hs::FLP}1, y[1] w[118]; P{ovoD1-18}3L P{FRTw[hs]}2A (for genotypes with wild-type Vinculin) or ΔVinc w[-]; P{hs::FLP}38/CyO; P{ovoD1-18}3L P{FRTw[hs]}2A (for genotypes with ΔVinc) males. Heat shocks were performed two times for 1 hr and 15 min each at 37°C at L1 and L2 larval stages. TalinIBS2-mCherry [31 (link)] was kindly provided by H.J. Bellen. The myosin heavy chain mutant used was Mhc[1] [45 (link)], kindly provided by S.I. Bernstein. IBS2-GFP recruitment to muscle attachment sites was performed with UAS::IBS2-GFP [15 (link)] expressed in muscles with P{Gal4-Mef2.R}3 (Bloomington Drosophila Stock Center).

To analyse the spectrum of mutation, we grouped the trios into higher taxonomic levels, that is, mammals, birds, fishes and reptiles. Thus, the percentages reported are based on the total candidate mutations from each group of species. We explored the genomic context of the mutations from a C or a G base to determine whether they were located in CpG sites (respectively followed by a G or preceded by a C) (see Supplementary Table 4). We phased the DNMs to their parental origin using the read-backed phasing method described previously (GitHub: https://github.com/besenbacher/POOHA)^{82 (link)}. This method uses the read-pairs containing a DNM and another heterozygous variant to determine the parental origin of the mutation when the heterozygous variant is present in both the offspring and one of the parents. The phasing allowed us to identify parental biases in the contribution of the DNMs by grouping multiple species to increase the number of phased mutations and obtain a minimum of 30 phased mutations per taxon. From this analysis, we omitted the Egyptian roussette (Rousettus aegyptiacus), Chinese tree shrew (Tupaia belangeri), griffon vulture (Gyps fulvus), blue-throated macaw (Ara glaucogularis), snowy owl (Bubo scandiacus) and Darwin’s rhea (Rhea pennata), as these could not be grouped with another monophyletic clade. To quantify the effect of parental age, a linear regression between the per-generation mutation rate and the average parental age at the time of reproduction was implemented using the lm function in R. Multiple linear regression was also used to identify whether paternal or maternal age was the strongest predictor of the empirical mutation rate.

Metagenomic DNA was isolated from approximately 180–200 mg of cecal digesta using genomic DNA columns (Macherey‐Nagel, Düren, Germany) according to Lagkouvardos et al. [38 (link)]. V3-V4 regions of the 16S rRNA genes were amplified using bacteria‐specific primers following a two-step procedure according to the Illumina sequencing protocol as described [38 (link)]. Amplicons were sequenced using a MiSeq system (Illumina, Inc., San Diego, CA, USA). Further processing of raw sequences was carried out as described recently [39 (link)]. Finally, sequences with a relative abundance > 0.1% in at least one sample were sorted, merged and operational taxonomic units (OTU) were picked at a threshold of 97% similarity. Taxonomic classification to the OTU was assigned using the SILVA database [40 (link)]. Further downstream analyses were done using Rhea (https://lagkouvardos.github.io/Rhea/). The differential abundance analysis of taxa was performed on the aggregated data at the different taxonomic levels as described [38 (link)]. For estimation of diversity within samples (α-diversity), the Shannon and Simpson indices, the most common indices to compare diversity, were calculated and transformed to the corresponding effective number of species according to Jost [41 (link)], because they are better suited at indicating the true diversity between samples and are minimally affected by the number of rare species. To measure the similarity between different microbial profiles, the β-diversity was determined by calculating generalized UniFrac distances with PERMANOVA statistical test as described previously [38 (link)]. Visualization of bacterial profiles among different groups was done by computation of non-metric multidimension distance scaling (NMDS) [42 (link)].