Introduction
Trematodes of the family Azygiidae Lühe, 1909 parasitise stomachs or body cavities in elasmobranchs and stomachs in freshwater teleosts and holosteans (Gibson Reference Gibson, Gibson, Jones and Bray2002). Representatives of the type genus of this family, Azygia Looss, 1899, possess a characteristic medium- or large-sized, elongated body, a small oral sucker, a large ventral sucker, and two large tandem testes. In the Russian Far East, these worms are known to infect mainly freshwater fishes, including such species as Esox reicherti Dybowski, 1869, Perccottus glenii Dybowski, 1877, Hucho taimen (Pallas 1773), and Channa argus (Cantor, 1842) (Mamaev, Oshmarin, Reference Mamaev, Oshmarin and Mamaev1971; Dvoryadkin Reference Dvoryadkin1977; Ermolenko et al. 1998; Besprozvannykh Reference Besprozvannykh2005; Vainutis et al. Reference Vainutis, Voronova, Mironovsky, Zhigileva and Zhokhov2023). The phylogenetic position of Azygiidae among other members of the Hemiurata is currently under debate. As a consequence, the taxonomic status of this family still remains unresolved. Taxonomists previously considered this group of trematodes as a separate suborder, Azygiata La Rue, Reference La Rue1957, or the order Azygiida Odening, 1963 (La Rue Reference La Rue1957; Skrjabin & Guschanskaya Reference Skrjabin, Guschanskaja and Skrjabin1958; Nagasawa et al. Reference Nagasawa, Urawa and Awakura1987; Littlewood Reference Littlewood2008; Sokolov & Zhukov Reference Sokolov and Zhukov2016). At present, most authors, based on the results of molecular phylogenetic analyses using ribosomal DNA gene sequence data (Olson et al. Reference Olson, Cribb, Tkach, Bray and Littlewood2003; Pérez-Ponce de León & Hernández-Mena Reference Pérez-Ponce de León and Hernández-Mena2019), recognize the status of this trematode group as a separate superfamily, Azygioidea Lühe, 1909 (Gibson Reference Gibson, Gibson, Jones and Bray2002; Olson et al. Reference Olson, Cribb, Tkach, Bray and Littlewood2003; Kostadinova & Pérez-del-Olmo Reference Kostadinova, Pérez-del-Olmo, Toledo and Fried2014; Pérez-Ponce de León & Hernández-Mena Reference Pérez-Ponce de León and Hernández-Mena2019), and as a member of the suborder Hemiurata Skrjabin & Guschanskaja, 1954. A phylogenetic analysis of Digenea based on complete mitochondrial sequence data and also using the whole mitochondrial DNA (mtDNA) genome of an azygiid representative, A. hwangtsiyui Tsin, 1933, obtained for the first time, has shown that Azygiidae represents a distinct branch, basal for most of trematode groups except Schistosomatidae Stiles & Hassall, 1898 (Wu et al. Reference Wu, Gao, Cheng, Xie, Yuan, Liu and Song2020). In our opinion, these results provide sufficient grounds for revising the taxonomic status of Azygiidae through further phylogenetic studies using mtDNA complete sequence data for different azygiid species. In our present study, we provide new data on the complete mtDNA sequence, inferred by the next-generation sequencing (NGS) approach, from two adult specimens of the trematode Azygia robusta Odhner, 1911 extracted from two salmonid species, the taimens Hucho taimen and Parahucho perryi, which were caught in two rivers of Primorsky Krai, Russia. This trematode species was earlier characterised morphologically by Besprozvannykh (Reference Besprozvannykh2005), who provided a detailed description of its life cycle. Our study aimed mainly to compare the structures and variations in the complete mitochondrial genomes of two azygiid species, analyse phylogenetic relationships using the new complete mtDNA sequence data on A. robusta, and interpret the obtained results to clarify the Azygiidae systematics.
Material and methods
Sample collection and DNA extraction
Adult worms were collected from the intestines of two naturally infected salmonids, a common taimen (H. taimen) caught in the Armu River (Besprozvannykh Reference Besprozvannykh2005) and a Sakhalin taimen (P. perryi (Brevoort, 1856)) caught in the Samarga River (unpublished, collected in 1987) (Table 1). The trematodes were killed with hot water and then fixed in 96% ethanol. Total DNA was extracted from the two specimens separately using a Qiamp Investigator kit (Qiagen, Germantown, MD, USA) according to the manufacturer’s protocol. Amount of total DNA was measured on a Qubit 3.0 fluorometer (Invitrogen, Waltham, MA, USA) and then used for NGS sequencing in a final amount of 100 ng.
Preparation of genome library for NGS
Libraries were prepared using an Ion Plus Fragment Library kit and unique adapters from an Ion Xpress Barcode Adaptors kit (ThermoFisher Scientific, Waltham, MA, USA) with pre-fragmentation on a Covaris M220 Focused-ultrasonicator (Covaris, LLC, Woburn, MA, USA). The preparation of polymerase chain reaction (PCR) emulsion and templates was done on an Ion One Touch 2 System (ThermoFisher Scientific) followed by sequencing on an Ion S5 sequencing platform using an Ion 540 chip.
The sequence quality and length distribution of raw reads were checked using FastQC 0.11.9 (Babraham Bioinformatics) and then the reads were assembled using SPAdes 3.14.1 (Nurk et al. Reference Nurk, Bankevich, Antipov, Gurevich, Korobeynikov, Lapidus, Prjibelski, Pyshkin, Sirotkin, Sirotkin, Stepanauskas, Clingenpeel, Woyke, JS, Lasken, Tesler, Alekseyev and Pevzner2013) with correction of IonTorrent data using the IonHammer tool available in the SPAdes software. The scaffolds containing mtDNA data were manually assembled in the MEGA X software (Kumar et al. Reference Kumar, Stecher, Li, Knyaz and Tamura2018).
Mitochondrial genome annotation was performed using the MITOS2 on-line software (Donath et al., Reference Donath, Jühling, Al-Arab, Bernhart, Reinhardt, Stadler, Middendorf and Bernt2009, available at http://mitos2.bioinf.uni-leipzig.de), and then the mitochondrial genome was manually assembled and aligned with that of A. hwangtsiyui (Wu et al. Reference Wu, Gao, Cheng, Xie, Yuan, Liu and Song2020) in MEGA X. Tandem repeats were searched using the Tandem Repeat Finders software (https://tandem.bu.edu/trf/trf.html). Search and analysis of the transfer RNA (tRNA) gene structure were performed in the ARWEN software (http://130.235.244.92/ARWEN/).
Codon usage, gene variations, and phylogenetic analyses
Alignments of nucleotide and amino acid sequences were performed by the ClustalW algorithm in MEGA X. Poorly aligned regions were removed using the Gblocks Server (http://phylogeny.lirmm.fr/phylo_cgi/one_task.cgi?task_type=gblocks).
Phylogenetic analysis was performed on the basis of concatenated amino acid sequences by the Maximum likelihood (ML) algorithm available in the PhyML 3.1 software (Guindon & Gascuel Reference Guindon and Gascuel2003) and by the Bayesian Inference (BI) method available in the MrBayes 3.2.6 software (Ronquist et al. Reference Ronquist, Teslenko, van der Mark, Ayres, Darling, Höhna, Larget, Liu, Suchard and Huelsenbeck2012). The ML algorithm was performed using an LG evolutionary model (Lee & Gascuel Reference Le and Gascuel2008), Subtree Pruning and Regrafting (SPR) tree topology search, and random sequence addition. The BI algorithm was performed using a protein model, a mixed set of substitution types, a mixed amino acid model, and uninformative amino acid substitution rates. The Monte Carlo Markov chains algorithm was performed with 1 000 000 generations during two independent runs, with sampling each 1000th generation and burning the first 25% of all generations. The average standard deviation of split frequencies was 0.000865, and that was enough for phylogenetic reconstruction. Significance of phylogenetic relationships was estimated with a posteriori probabilities (Huelsenbeck et al. Reference Huelsenbeck, Ronquist, Nielsen and Bollback2001) for the BI algorithm and an approximate likelihood-ratio test (Anisimova & Gascuel Reference Anisimova and Gascuel2006) for the ML algorithm. Codon usage statistics was calculated for concatenated protein-coding gene sequence data in MEGA X. Analysis of correlation between the number of variable sites and the gene length was performed using Pearson’s correlation coefficient in Statistica 13 software (TIBCO Software Inc. 2017).
Phylogenetic relationships were inferred using sequences of our samples and other trematode species accessed from the NCBI GenBank database (Table 1). The two annotated mitochondrial genomes have been deposited in GenBank under accession numbers OR350239 and OR350240, while raw Sequence Read Archive (SRA) sequencing data are available under accession numbers SAMN36469092–SAMN36469093.
Results
Brief visual morphological identification of the species
In this study we first performed a brief visual morphological analysis of the paratypes of trematodes used for the NGS analysis. The general view of the azygiid worms from H. taimen caught in the Armu River and from P. perryi caught in the Samarga River can be seen in Figures 1 and 2, respectively. Both specimens possess the main diagnostic characteristics of A. robusta, including the round pharynx and vitellaria extending beyond the posterior end of the second testis to half the distance between the second testis and the posterior end of the body (Skrjabin & Guschanskaja Reference Skrjabin, Guschanskaja and Skrjabin1958; Bauer Reference Bauer1987). These morphological characteristics were observed clearly in both specimens (Figures 5, 6). For the NGS analysis, we used the mature trematode specimens that had been found in H. taimen from the Armu River cultivated from cercariae, identified as A. robusta, and published in the study of Besprozvannykh (Reference Besprozvannykh2005), who described the life cycle of this trematode species. Thus, we keep the name Azygia robusta for both trematode specimens used for the NGS analysis.
Sequence quality and coverage
We obtained 3.5–4.5 million reads for two specimens of A. robusta. The sequence quality after FastQC was acceptable. Phred 33 values were 20–30 (mode 26) and decreased slightly in long reads. The sequence length was 25–449 bp; for most reads, the length was 120–240 bp. The GC content and numbers of duplications and adapters did not exceed the norm. The mean coverage across the mitochondrial DNA was 103X and 204X for the two specimens of A. robusta.
General characteristics of the Azygia robusta mitochondrial genome
The mitochondrial genome of A. robusta had a length of 13 857 bp and contained 12 protein-coding genes, two ribosomal genes, 22 tRNA genes, and two non-coding regions: short (SNCR) and long (LNCR) (Figure 3, Table 2). Alternative read variants were absent from the NGS raw data, and no intraspecific variable positions were observed. The nucleotide composition in the A. robusta mitochondrial genome was as follows: A, 16.5%; T (U), 40.9%; C, 14.4%; and G, 28.2%. The nucleotide pair frequency was 57.4% for the AT-content and 42.6% for the GC-content, showing a bias towards T over A (AT skew = -0.43) and G over C (CG skew = 0.32), respectively.
Protein-coding genes of mitochondrial genome
The total length of the sequences of 12 protein-coding genes in the complete mitochondrial genome was 10 110 bp. The arrangement of protein-coding genes was as follows: cox3–cytb–nad4L–nad4–atp6–nad2–nad1–nad3–cox1–cox2–nad6–nad5. The start-codons for protein-coding genes were ATG or GTG, except the cox1 gene that started with TTG codon, as well as those for A. hwangtsiyui and nad3 gene that started with GGT codon (Table 2). The nucleotide composition of the assembled protein-coding part of the mitochondrial genome sequence was as follows: A, 14.5%; T (U), 43.3%; C, 14.0%; and G, 28.2%. The nucleotide pair frequency was 57.8% for the AT-content and 42.2% for the GC-content, showing a bias towards T over A (AT skew = -0.5) and G over C (CG skew = 0.34).
LNCR, long non-coding region; SNCR, short non-coding region.
* tRNA missed the DHU-arm
** tRNA missed the T-stem
Codon usage statistics for A. robusta were consistent with the proportions in the nucleotide composition: the most common triplets contained T (U) and/or G bases, namely UUU (with a frequency of 9.79%), UUG (5.99%), GUU (5.88%), GGG (4.4%), UGU (3.82%), and GUG (3.39%).
A total of 3 386 amino acids were encoded by the mitochondrial protein-coding genes in A. robusta. Of these, a maximal frequency was observed for leucine (15.0%), valine (12.2%), and phenylalanine (11.7%); the frequencies for lysine (1.18%) and glutamine (1.21%) were minimal compared to other amino acids. The amino acid frequencies of the mitochondrial protein sequences of A. hwangtsiyui were similar to those of A. robusta, with no marked differences observed (Table 3).
Intra- and interspecific variation of complete mtDNA sequences
Overall, the nucleotide sequences of the complete mitochondrial genomes, including all genes and non-coding fragments, of the two A. robusta specimens differed from each other by 0.12 ± 0.03%. The concatenated protein-coding gene sequences between these two specimens differed by 0.11 ± 0.03%, and amino acid sequences by 0.06 ± 0.05%. Six of 12 protein-coding genes demonstrated intraspecific variation in A. robusta (Table 4); a total of 12 substitutions were revealed, with each gene containing from one to three variable sites, transitions T/C (67%) or A/G (25%), and a single transversion T/G in nad6 gene.
The difference between the nucleotide sequences of the complete mitochondrial genomes of A. robusta and A. hwangtsiyui was 26.95 ± 0.35%; between the concatenated protein-coding nucleotide sequences, 26.00 ± 0.43%; and between the amino acid sequences, 30.15 ± 0.82%. The interspecific variation of protein-coding genes between A. robusta and A. hwangtsiyui ranged from 20.5 ± 0.9% (cox1) to 30.7 ± 1.2% (nad5) (Table 4). The results of correlation analysis using Pearson’s correlation coefficient showed a high positive correlation (r = 0.95) between the number of variable sites and the gene length for interspecific comparison of protein-coding gene variations in the two Azygia species (Figure 4).
Phylogenetic analysis
The maximum likelihood (ML) and Bayesian Inference (BI) algorithms were based on alignment of 2 280 amino acids available after Gblocks processing. Overall, mitochondrial genomes of 62 species, including 61 digenean and one cestode species, Diphpyllobothrium latum (Linnaeus, 1758) Lühe, 1899, were incorporated into the phylogenetic analysis. As the BI tree topology showed, the digeneans could be subdivided into three highly supported clades (Figure 5). The first clade was early divergent and consisted of three Azygia specimens: one A. hwangtsiyui (GenBank accession no. MN844889) and two A. robusta (our material). The second clade represented the order Diplostomida, including species of the families Schistosomatidae Stiles & Hassal, 1898, Clinostomidae Lühe, 1901, Cyathocotylidae Mühling, 1898, and Brachylaimidae Joyeux & Foley, 1930. The third clade comprised 47 species from 18 families, representing seven suborders. The topology of this subclade completely agreed with the previous phylogenetic reconstructions of Digenea (Ivashko et al. Reference Ivashko, Semenchenko, Solodovnik and Atopkin2022). The suborder Xiphidiata Olson, Cribb, Tkach, Bray & Littlewood, Reference Olson, Cribb, Tkach, Bray and Littlewood2003 was polyphyletic and appeared as two independent groups. The first group consisted of Brachycladium goliath (van Beneden, 1858) and species of the genus Paragonimus Braun, 1988 (Xiphidiata). This group was closely related to Opisthorchiata. Plagiorchis maculosus (Rudolphi, 1802) appeared as sister to the above-mentioned Xiphidiata and Opisthorchiata species. The second group of Xiphidiata included representatives of Dicrocoeliidae Odhner, 1911 (Dicrocoelium dendriticum (Rudolphi, 1819), D. chinensis (Sudarikov and Ryjikov, 1951) Tang and Tang, 1978, and Eurytrema pancreaticum (Janson, 1899)), Eucotylidae Skrjabin, 1924 (Tamerlania zarudnyi Skrjabin, 1924), and Prosthogonimidae Lühe, 1909 (Prosthogonimus cuneatus (Rudolphi, 1802) Braun, 1901). This group appeared as sister to the subclade that contained species of the suborder Pronocephalata Olson, Cribb, Tkach, Bray, Littlewood, Reference Olson, Cribb, Tkach, Bray and Littlewood2003. The suborder Haploporata Pérez-Ponce de León & Hernández-Mena, Reference Pérez-Ponce de León and Hernández-Mena2019, represented by Parasaccocoelium mugili Zhukov, 1971 and Carassotrema koreanum Park, 1938, formed a basal branch within the third clade.
In general, the ML tree topology was similar to that of the BI tree, demonstrating three main clades within Digenea (Figure 6). The marked differences from the BI tree topology were as follows: (1) Dicrocoeliidae (Xiphidiata) formed a separate basal subclade within the third clade, creating polyphyly for Xiphidiata, and (2) representatives of Haploporata, C. koreanum and P. mugili, were within a single subclade with Plagiorchis maculosus, with this subclade being sister to closely related representatives of Xiphidiata (Paragonimus spp. and Brachycladium goliath) and Opisthorchiata. These differences between the BI and ML tree topologies were also reported in previous studies (Atopkin et al. Reference Atopkin, Semenchenko, Solodovnik, Ivashko and Vinnikov2021; Ivashko et al. Reference Ivashko, Semenchenko, Solodovnik and Atopkin2022).
Discussion
Mitochondrial genome variations in A. robusta and A. hwangtsiyui
The complete mitochondrial genome structure of A. robusta was highly similar to that of A. hwangtsiyui in gene arrangement, the existence of two non-coding regions separated from each other by tRNA-Gly (G) gene, and a higher level of AT content relative to GC content for both mitochondrial genome sequences and protein-coding genes. Also, the two azygiid species had the same most frequent codons and the same start-codon (TTG) for the cox1 gene. A difference was revealed in the start-codon of the nad3 gene, which started with GGT in A. robusta vs. ATG in A. hwangtsiyui. There were also differences in the lengths of some protein-coding genes between A. robusta and A. hwangtsiyui: 1275 vs. 1272 bp, respectively, for the nad4 gene; 504 vs. 513 bp for atp6; 903 vs. 906 bp for nad1; 1548 vs. 1564 bp for cox1; 450 vs. 444 bp for nad6; and 1598 vs. 1600 bp for nad5.
New definitive host of Azygia robusta from the Russian Far East
To date, four species of definitive hosts for trematodes Azygia spp. are known from the Russian Far East: the northern snakehead Channa argus (Cantor, 1842) and the Amur pike Esox reicherti Dybowski, which are freshwater fishes, and the common taimen Hucho taimen (Pallas, 1773) and the Chinese sleeper Perccottus glenii Dybowski, 1877, which are freshwater/brackish-water fishes (Mamaev, Oshmarin, Reference Mamaev, Oshmarin and Mamaev1971; Dvoryadkin Reference Dvoryadkin1977; Ermolenko et al. 1998; Besprozvannykh Reference Besprozvannykh2005; Vainutis et al. Reference Vainutis, Voronova, Mironovsky, Zhigileva and Zhokhov2023). In this study, we have extended the list of definitive hosts for this region by adding the Sakhalin taimen, Parahucho perryi (Brevoort, 1856). This fish is one of the world’s largest salmonids, reaching maturity at 6–8 years of age and living for more than 20 years. The species’ geographic range is confined to the Sea of Japan, from the southern Kuril Islands and Primorsky Krai, Russia, to Hokkaido, Japan. Parahucho perryi occupy a variety of habitats including upper and lower reaches of rivers, lakes, brackish-water lagoons, estuaries, and coastal marine waters (Zolotukhin & Semenchenko Reference Zolotukhin and Semenchenko2008; Fukushima et al. Reference Fukushima, Shimazaki, Rand and Kaeriyama2011). The ecological features of P. perryi are favorable for the trematode A. robusta to infect this fish species. Moreover, one of the azygiid species, A. perryi, has been reported as a parasite for P. perryi from Japan (Nagasawa et al. Reference Nagasawa, Urawa and Awakura1987; Popiołek et al. Reference Popiolek, Kusznierz, Kotusz and Witkowski2013), while A. robusta is known to parasitise salmonids (Nikolić et al. Reference Nicolić, Bilbija, Nedic, Simonovic and Djikanovic2018). Thus, we here provide the first record of a new definitive host, Parahucho perryi, for the trematode Azygia robusta from the Russian Far East.
Systematics and phylogenetic relationships of Azygiidae
The systematic position of Azygiidae is still unclear because of the lack of molecular data for representatives of this group, and, as a consequence, controversies arise between interpretations of morphological and molecular data. Most authors recognize the status of this trematode group as a separate superfamily (Gibson Reference Gibson, Gibson, Jones and Bray2002; Olson et al. Reference Olson, Cribb, Tkach, Bray and Littlewood2003; Kostadinova & Pérez-del-Olmo Reference Kostadinova, Pérez-del-Olmo, Toledo and Fried2014; Pérez-Ponce de León & Hernández-Mena Reference Pérez-Ponce de León and Hernández-Mena2019). However, viewpoints on the membership of Azygioidea at a higher taxonomic level are different. These worms were considered as a separate suborder, Azygiata (La Rue Reference La Rue1957; Skrjabin & Guschanskaya Reference Skrjabin, Guschanskaja and Skrjabin1958), or the order Azygiida Odening, 1963 (Littlewood Reference Littlewood2008; Sokolov & Zhukov Reference Sokolov and Zhukov2016). At present, this superfamily is recognized as a member of the suborder Hemiurata mainly on the basis of data inferred from molecular phylogenetic analyses using ribosomal DNA gene sequences (Olson et al. Reference Olson, Cribb, Tkach, Bray and Littlewood2003; Pérez-Ponce de León & Hernández-Mena Reference Pérez-Ponce de León and Hernández-Mena2019). However, as the latest studies have shown, the Azygiida is a valid order (Ramilo et al. Reference Ramilo, Abollo and Pascual2023). The first complete mitochondrial genome sequences for a representative of Azygiidae, Azygia hwangtsiyui, were obtained by Wu et al. (Reference Wu, Gao, Cheng, Xie, Yuan, Liu and Song2020). These data were applied to the reconstruction of the phylogenetic position of Azygiidae within Digenea using a dataset of concatenated amino acid sequences representing 12 protein-coding genes. The position of A. hwangtsiyui (Azygiidae) was considered the ‘most basal lineage of the Digenea’; however, in that study, Schistosomatidae rather than Azygiida was basal for Digenea (Wu et al. Reference Wu, Gao, Cheng, Xie, Yuan, Liu and Song2020). Nevertheless, the authors did not provide any definitive conclusion about the systematic position of Azygiidae and stated that ‘the family Azygiidae still awaits investigation of relationships based on a much wider taxon sampling and more mitogenome datasets’ (Wu et al. Reference Wu, Gao, Cheng, Xie, Yuan, Liu and Song2020). Our results clearly demonstrate that this statement is relevant. The introduction of one new azygiid species into the phylogenetic analysis based on concatenated amino acid sequences of 12 protein-coding mitochondrial genes has considerably changed the phylogenetic position of Azygiidae within Digenea. In contrast to the results from previous studies, the present phylogenetic tree consists of three main, highly supported digenean clades, including the early diverging clade Azygiidae and two sister clades, Diplostomida and other digeneans. On the one hand, this result supports the taxonomic status of Azygiidae as a separate order within Digenea, which largely agrees with the previous results by Wu et al. (Reference Wu, Gao, Cheng, Xie, Yuan, Liu and Song2020) that showed a basal position of Azygiidae relative to other Digenea, except Schistosomatidae. On the other hand, our results do not confirm the hypothesis, advanced in our previous studies, about the consistency between phylogenetic relationships and gene rearrangement within mitochondrial genomes of Schistosomatidae and other trematodes (Atopkin et al. Reference Atopkin, Semenchenko, Solodovnik, Ivashko and Vinnikov2021; Ivashko et al. Reference Ivashko, Semenchenko, Solodovnik and Atopkin2022). In particular, the basal position of Azygiidae relative to other trematodes and the emergence of Cyathocotyle prussica and Clinostomum complanatum within a single clade with Schistosomatidae are evidence against this hypothesis. However, in this respect, we agree with Wu et al. (Reference Wu, Gao, Cheng, Xie, Yuan, Liu and Song2020), who indicated the need for additional data, with complete mitochondrial genome sequences obtained for more Azygiidae species and other unstudied digenean taxa, to provide a sufficient basis for conclusions about the systematics of this family.
Acknowledgements
We are deeply thankful to Vladimir V. Besprozvannykh, Dr. Sci. Biol., the head of the Department of Parasitology, Federal Scientific Center of East Asian Terrestrial Biodiversity FEB RAS, for providing the material on Azygia robusta (the specimens for the DNA extraction and the microscope slides of paratypes for the morphological analysis) and to Dr. Kirill A. Vinnikov, the Director of the Institute of the World Ocean, Far Eastern Federal University, for providing the laboratory equipment and constructive comments while performing NGS.
Financial support
This study was supported by the Russian Scientific Foundation, project number 22-24-00896.
Competing interest
None.
Ethical standard
All applicable institutional, national and internationalguidelines for the care and use of animals were followed.