Abstract

Simple Summary

   Interspecific hybrids of F1 generations have frequently been
   characterized by high vigour resulting from heterosis advantage. In
   contrast, post-F1 generations are expected to express hybrid breakdown,
   i.e., they suffer from low viability and survival, reproductive
   abnormalities or sterility and limited ecological performance.
   Resistance or susceptibility to parasites is one of the measures
   reflecting hybrid vigour. The present study aimed to analyse the
   experimental infection of the blood-feeding generalist parasite
   Paradiplozoon homoion (Monogenea) in two target fish species, Abramis
   brama and Rutilus rutilus, and their reciprocal F1 hybrids and
   backcross hybrids, and to reveal potential parasite-induced changes in
   their transcriptome profiles of head kidney. We hypothesized various
   effects of hybridization in terms of parasitism in F1 hybrids and
   backcross hybrids reflected by differential gene expression. The number
   of differentially expressed genes (DEGs) differed between fish lines
   with a lower number of DEGs in F1 hybrids and a higher number in
   backcross hybrids when compared to the parental lines, A. brama and R.
   rutilus. Backcross hybrids were more infected than F1 hybrids and
   parental lines. DEG analyses revealed the role of heme binding,
   erythrocyte differentiation and immunity-related processes in fish
   after infection by blood-feeding P. homoion. Using GO and KEGG
   analyses, we revealed the similarity in DEGs between two backcross
   generations of hybrids. This finding may indicate a potential
   consequence of hybrid breakdown in backcross generations. Gene
   expression in less parasitized F1 hybrids is in line of hybrid
   advantage.

Abstract

   Hybrid generations usually face either a heterosis advantage or a
   breakdown, that can be expressed by the level of parasite infection in
   hybrid hosts. Hybrids are less infected by parasites than parental
   species (especially F1 generations) or more infected than parental
   species (especially post-F1 generations). We performed the experiment
   with blood-feeding gill parasite Paradiplozoon homoion (Monogenea)
   infecting leuciscid species, Abramis brama and Rutilus rutilus, their
   F1 generation and two backcross generations. Backcross generations
   tended to be more parasitized than parental lines and the F1
   generation. The number of differentially expressed genes (DEGs) was
   lower in F1 hybrids and higher in backcross hybrids when compared to
   each of the parental lines. The main groups of DEGs were shared among
   lines; however, A. brama and R. rutilus differed in some of the top
   gene ontology (GO) terms. DEG analyses revealed the role of heme
   binding and erythrocyte differentiation after infection by
   blood-feeding P. homoion. Two backcross generations shared some of the
   top GO terms, representing mostly downregulated genes associated with
   P. homoion infection. KEGG analysis revealed the importance of
   disease-associated pathways; the majority of them were shared by two
   backcross generations. Our study revealed the most pronounced DEGs
   associated with blood-feeding monogeneans in backcross hybrids,
   potentially (but not exclusively) explainable by hybrid breakdown. The
   lower DEGs reported in F1 hybrids being less parasitized than backcross
   hybrids is in line with the hybrid advantage.

   Keywords: Monogenea, Paradiplozoon homoion, freshwater fish,
   hybridization, RNA seq, differential gene expression, hybrid heterosis,
   hybrid breakdown

1. Introduction

   Monogeneans are mainly fish ectoparasites, frequently infecting gills,
   fins and skin and exhibiting a direct life cycle (i.e., an intermediate
   host is not involved in their life cycle). Their short generation time
   can result in exponential population growth. Serious effects on the
   health condition of wild-living fish have rarely been found; however,
   the considerable economic impact of monogenean infection on the health
   of aqua-cultured or farmed fish species has been well documented, i.e.,
   monogeneans reduce fish growth and may cause fish morbidity and
   mortality [[34]1,[35]2,[36]3,[37]4,[38]5,[39]6].

   Paradiplozoon (Diplozoidae) are relatively large-bodied monogenean
   parasites (with a body size ranging in different species from 0.8 to
   8.2 mm) and common gill parasites on cyprinoid fish [[40]7]. Among
   Paradiplozoon species, Paradiplozoon homoion is a generalist species
   with the highest range of host species [[41]8,[42]9,[43]10,[44]11].
   Paradiplozoon homoion, like other diplozoid monogeneans, is
   representative of blood-feeding ectoparasites that degrade host blood
   components predominantly intracellularly [[45]12]. Usually, the
   intensity of infection by P. homoion (similarly to other diplozoid
   species) is low and its prevalence varies depending on the reproduction
   period of the parasite, the localities of sampling and the host species
   (e.g., [[46]13,[47]14,[48]15,[49]16]). Adult diplozoids are obligatory
   blood-feeders, and some diplozoid species, e.g., Eudiplozoon
   nipponicum, are known as important pathogens of farmed cyprinoid fishes
   [[50]7,[51]8,[52]17,[53]18].

   Two fish species, common bream (Abramis brama) and roach (Rutilus
   rutilis), both representatives of the family Leuciscidae, are
   widespread in European freshwaters. The hybridization of these species
   was previously widely documented [[54]13,[55]19,[56]20,[57]21,[58]22].
   Hybrids of A. brama and R. rutilus in natural habitats primarily result
   from reciprocal crosses of parental species (i.e., hybrids of F1
   generations are present), whilst the presence of post-F1 hybrids is
   negligible [[59]20,[60]21] indicating that F1 hybrids have reproductive
   disadvantages, or post-F1 hybrids (backcross or F2) experience
   breakdown, or at least low survival. Both A. brama and R. rutilus as
   well as their F1 hybrids are susceptible to Paradiplozoon homoion
   [[61]13], although R. rutilus is documented as the most common host
   species for P. homoion
   ([62]https://www.nhm.ac.uk/research-curation/scientific-resources/taxon
   omy-systematics/host-parasites/database, accessed on 30 May 2023).

   Generally, hybrids of F1 generations have frequently been characterized
   by high hybrid vigour, i.e., heterosis advantage, which arises when
   traits are additive, or due to overdominance resulting in synergic
   effect of parental alleles on hybrid vigour. F1 hybrids often exhibit
   superiority over their parents, e.g., greater growth, survival and
   environmental tolerance or resistance to parasites (for fish that have
   been documented, see, e.g., [[63]13,[64]23,[65]24,[66]25,[67]26]). For
   the F1 generation of evolutionarily divergent A. brama and R. rutilus,
   hybrid heterosis was revealed through the utilization of a broader
   trophic spectrum and a higher tolerance to fluctuations in food supply
   compared to parental species [[68]20,[69]27], and in terms of greater
   survival at early developmental stages [[70]28], faster growth [[71]29]
   and lower susceptibility to parasites when compared to parental species
   [[72]13,[73]16,[74]30]. Studies investigating the metazoan parasite
   communities in A. brama, R. rutilus and their F1 hybrids either in
   natural habitats or under experimental conditions consistently reported
   lower monogenean abundance in F1 hybrids when compared to parental
   species, which is in line with hybrid heterosis generally predicted for
   F1 generations [[75]13,[76]16]. Šimková et al. [[77]30] showed that the
   diversity of MHC II (major histocompatibility complex) genes, i.e.,
   highly polymorphic genes involved in the adaptive immunity of
   vertebrates, in F1 hybrids of A. brama and R. rutilus was intermediate
   between two leuciscid species, supporting hybrid advantage with respect
   to coping with parasites (i.e., hybrids with intermediate number of MHC
   alleles and intermediate positively selected sites potentially
   reflecting optimal MHC carried fewer parasites than either parental
   species). However, post-F1 generations of hybrids are expected to
   express so-called hybrid breakdown, i.e., their biological fitness is
   reduced as they suffer from low viability, low survival, reproductive
   abnormalities or sterility and limited ecological performance. The
   phenomenon of hybrid breakdown results from various genetic
   incompatibilities [[78]31,[79]32]; hybrids exhibit many disadvantageous
   traits, often resulting from the disruption of gene expression
   regulation or the disruption of nuclear–mitochondrial gene interactions
   [[80]33]. Hybrid breakdown in fish was revealed by several studies
   (e.g., [[81]25,[82]34,[83]35,[84]36,[85]37]). Regarding A. brama and R.
   rutilus hybridization, Stolbunova et al. [[86]38] showed the
   differences in morphology and suggested the differences in viability
   (based on asymmetrical ratio of offspring genotypes) between
   backcrossed hybrids with the mtDNA of R. rutilus and mtDNA of A. brama,
   which was interpreted as evidence of varying degrees of cyto-nuclear
   compatibility of the genomes of A. brama and R. rutilus.

   However, the experimental study of monogenean parasite infection in A.
   brama, R. rutilus and their various hybrid lines (including F1 and
   post-F1 generations) did not support hybrid breakdown (the expectation
   of high parasite abundance in post-F1 generations) and, in contrast,
   showed that according to the abundance of host-specific parasites,
   backcross hybrid generations exhibited similarities with parental
   species whose genes contributed more to the backcross genotype
   [[87]16].

   Knowledge about the immune response in fish hosts to monogenean
   infection in terms of changes in gene expression is still
   underexplored. However, several recent studies have focused on the
   immune response of fish to monogenean infection, especially to
   Gyrodactylus infection (viviparous Monogenea) and Dactylogyrus
   infection (oviparous Monogenea), underlining the role of mucosal
   immunity [[88]39], the response through immune gene expression in
   spleen or head kidney [[89]40,[90]41], or in both skin and/or gills
   following infection by Gyrodactylus [[91]42,[92]43], or the coinfection
   of viviparous and oviparous representatives of monogeneans, i.e.,
   Gyrodactylus and Cichlidogyrus [[93]44]. In general, there are not a
   lot of studies investigating hybrid expression profiles in relation to
   parasitism, especially such studied are limited in fish, which are,
   however, organisms of high economic interest. Until now, no studies
   focusing on the effects of monogenean parasite infection on
   differential gene expression in the immune organs of pure species and
   their intergeneric hybrids have been performed, and, at the same time,
   there are no studies documenting potential changes in immune gene
   expression after diplozoid infection.

   In the present study, we focused on the common generalist monogenean
   parasite P. homoion, the presence of which was previously demonstrated
   in two target fish species, A. brama and R. rutilus and their F1
   hybrids in natural habitats [[94]13] and in their backcross hybrids
   under experimental infection [[95]16]. We performed an experimental
   study to compare the susceptibility of species lines (A. brama and R.
   rutilus), F1 hybrids and backcross hybrids using the same infectious
   dose of the P. homoion larval stage (free-living oncomiracidium); we
   also investigated transcriptome profile changes between infected and
   control fish of pure breed and hybrid lines. We hypothesized various
   effects of hybridization in terms of parasitism in F1 hybrids and
   backcross hybrids reflected by differential gene expression.
   Specifically, we expected differences in parasite load and
   transcriptome profile response between F1 hybrids and backcross
   hybrids, hypothesizing hybrid advantage for F1 hybrids and hybrid
   breakdown for backcross hybrids.

2. Material and Methods

2.1. Experimental Fish Lines

   Specimens of A. brama, R. rutilus and their natural F1 hybrids with A.
   brama mtDNA were collected from the Hamry Reservoir (49.73724 N,
   15.91395 E; the Czech Republic). The identification of hybrids was
   based on meristic traits (i.e., the number of gill rakers, the number
   of scales in the lateral line and the number of branched rays in the
   anal fin) and molecular markers (the cytochrome b gene and
   microsatellite loci following Krasnovyd et al. [[96]13]). Artificial
   breeding of the fish transported to the breeding facility was performed
   as described in Dedic et al. [[97]16]. Hormonal stimulation of fish for
   ovulation/spermiation followed Gela et al. [[98]45]; for details, see
   also Dedic et al. [[99]16]. Oocytes of ovulating females and sperm were
   sampled according to Gela et al. [[100]45] and Linhart et al.
   [[101]46]. The artificial breeding was performed for three pairs in
   each parental combination (see below); the offspring of different
   parental combinations representing the same fish line was mixed and 5
   specimens per each fish line were randomly selected for the experiment.

   The following parental combinations were used for artificial breeding
   ([102]Figure 1): (1) A. brama female and A. brama male, (2) R. rutilus
   female and R. rutilus male, (3) A. brama female and R. rutilus male,
   (4) F1 hybrid female (with A. brama mtDNA) and A. brama male and (5) F1
   hybrid female (with A. brama mtDNA) and R. rutilus male. Because of the
   limited number of free-living larval stage, oncomiracidium, used for
   experimental infection (see below), only F1 generation of hybrids with
   A. brama mtDNA, which is also naturally present in freshwaters, was
   used in this study. Two years old fish of five lines corresponding to
   R. rutilus line, A. brama line, F1 hybrid line with A. brama mtDNA
   (termed F1 A. brama × R. rutilus) and two backcross lines, i.e., the
   backcross generation resulting from the crosses of both parents with A.
   brama mtDNA (F1 hybrid female with A. brama mtDNA and A. brama male,
   termed backcross of F1 hybrid × A. brama) and the backcross generation
   resulting from the crosses of parents with different mtDNA (F1 hybrid
   female with A. brama mtDNA and R. rutilus male, termed backcross of F1
   hybrid × R. rutilus) were used for the experiment. Standard fish length
   in millimetres (shown as mean ± SD) for each fish line was measured: R.
   rutilus 79.58 ± 3.68, A. brama 82.67 ± 2.99, F1 A. brama × R. rutilus
   81.42 ± 7.96, backcross of F1 hybrid × R. rutilus 80.5 ± 9.84 and
   backcross of F1 hybrid × A. brama 91.17 ± 9.36.

Figure 1.

   [103]Figure 1
   [104]Open in a new tab

   Schema of artificial breeding. Five fish lines used for the experiment
   are shown (A. brama, R. rutilus, F1 A. brama × R. rutilus, backcross of
   maternal F1 hybrid with A. brama mtDNA and paternal A. brama, backcross
   of maternal F1 hybrid with A. brama mtDNA and paternal R. rutilus).
   Pure species lines are shown in blue, F1 hybrids are shown in orange,
   backcross hybrids are shown in green.

2.2. Monogenean Infection

   A total of 5 biological replicates (i.e., 5 fish specimens) per fish
   line were used in the experiment. Each naïve fish specimen was placed
   individually in a small aquarium (10 L) and acclimatized for one week
   to experimental room conditions. Temperature of the room was adjusted
   to 20 °C and aeration in each small aquarium was assured by an air
   stone roller of size 1 cm.

   Eggs of P. homoion were collected from the gills of donor fish
   specimens (donor Rutilus rutilus reared under aquarium conditions) and
   10–20 eggs were placed into embryo dishes (type 546/40 × 40, VITRUM
   Rožnov, Czech Republic) filled with tap water that was allowed to stand
   for 24 h prior to egg transfer. The embryo dishes were then placed into
   a laboratory incubator at 25 °C. The condition of the eggs was
   monitored twice daily using a stereomicroscope (Olympus SZX 7, Tokyo,
   Japan); the water was also changed at the same time. Hatching occurred
   after 2–3 days. The hatched oncomiracidia were used to infect naïve
   fish specimens. Twenty-five oncomiracidia specimens were transferred
   into the aerated and standing tap water in each aquarium in which a
   single naïve fish specimen was placed. Aeration was stopped for 2 h
   immediately after oncomiracidia were released into the water to ensure
   that the oncomiracidia reached the fish (following Pečínková et al.
   [[105]47], who reported the presence of successful oncomiracidia on
   fish 2 h from their introduction). Subsequently, aeration was run to
   ensure fish survival during the experiment. Fish were investigated for
   the presence of parasites one month after oncomiracidia infection
   following Pečínková et al. [[106]47] who showed that the development of
   P. homoion from egg to sexually mature adult stage took approximately
   33 days. All fish were dissected, i.e., fish were gently sacrificed
   following law n. 246/1992 of the Czech Republic (Act of the Czech
   National Council for the Protection of Animals Against Cruelty); then,
   the gills were removed and using stereomicroscope checked for the
   presence of P. homoion and head kidney of fish was removed using
   sterile instruments. All parasite specimens of P. homoion were counted.
   Head kidney was collected by fish dissection from individual fish
   including 5 non-infected control specimens per each fish line and 5 P.
   homoion-infected specimens per fish line. Tissue of each fish was
   separately submerged in Ambion RNAlater stabilization solution (Thermo
   Fisher Scientific, Waltham, MA, USA). The tubes with tissues were
   stored at −80 °C until the isolation of total RNA.

2.3. RNA Extraction and Library Preparation

   Total RNA was isolated from the head kidney of each fish specimen (25
   in total). For extraction, PureLink RNA Mini Kit (Ambion) with Trizol
   reagent (Thermo Fisher Scientific) and on-column PureLink DNase
   treatment were used according to the manufacturer’s protocol. Reagent
   and buffer volumes were adjusted according to the weight of tissue
   entering the isolation process (15.3 mg on average). The final elution
   was performed using 100 µL of RNAse-free water in the first step and
   the primal eluate in the second step. The yield and concentration of
   RNA isolates were checked using a Qubit^TM 4 fluorometer (Invitrogen by
   Thermo Fisher Scientific) and Qubit RNA HS Assay Kit (Thermo Fisher
   Scientific). The quality and integrity of RNA were analysed using RNA
   6000 Nano Kit on a 2100 Bioanalyzer instrument (Agilent Technologies,
   Santa Clara, CA, USA).

   All RNA isolates were normalized by dilution at a uniform concentration
   of 20 ng/µL with RNase-free water. They served as templates for DNA
   library preparation and for the reverse transcription of total RNA into
   single-stranded cDNA using High Capacity RNA-to-cDNA Kit (Applied
   Biosystems by Thermo Fisher Scientific) in twice the reaction volume
   recommended by the manufacturer.

   All 25 fish samples (RNA integrity number—RIN > 7) were used for DNA
   library preparation. We used 500 ng of total RNA for mRNA enrichment
   using the Poly(A) mRNA Magnetic Isolation Module (New England Biolabs).
   Subsequently, NEBNext^® Ultra™ Directional RNA Library Prep Kit for
   Illumina^® and NEBNext^® Multiplex Oligos for Illumina^® (Dual Index
   Primers Set 2, New England Biolabs) were used for library preparation,
   with 11 PCR cycles utilized for PCR enrichment. RNA fragmentation (13
   min at 94 °C) and the size selection conditions (a bead volume of 30 µL
   and 15 µL for the first and second bead selections, respectively) were
   further modified in the protocol. The quantification of DNA libraries
   was performed on a Qubit^TM 4 fluorometer (Invitrogen by Thermo Fisher
   Scientific) using Qubit dsDNA HS Assay Kit, and quality and size
   control were performed on a 2100 Bioanalyzer with DNA 1000 Kit (Agilent
   Technologies). Finally, amplicons were pooled in equimolar amounts. The
   final concentration of each particular library in the pool was 10 nM in
   the pool. Libraries were sequenced by Macrogen Europe B.V. (Amsterdam,
   The Netherlands) on an Illumina NovaSeq6000 system (one line on an S4
   flow cell) in a paired-end configuration producing 150 bp long reads.

2.4. NGS Data Analyses

   A quality check of raw paired-end fastq reads was carried out by FastQC
   [[107]48] and a contamination check against human, mouse, yeast, E.
   coli and other organisms by BioBloom tools [[108]49]. The quality and
   Illumina adapter trimming of raw reads was performed using Trimmomatic
   [[109]50]. As a preliminary step, all samples were aligned to the
   Carassius auratus (ASM336829v1-104) and Cyprinus carpio (common carp;
   GCA_000951615.2) genomes using STAR [[110]51]. The mapped reads were
   deduplicated by Picard’s MarkDuplicates [[111]52]. The quality of
   mapping was checked using Samtools [[112]53] for overall statistics and
   Picard for normalized gene coverage, aligned reads strandedness
   correctness and aligned reads assignment in the reference. The
   quantification of gene expression was performed by featureCounts
   [[113]54]. All compatible results and statistics were processed by
   MultiQC [[114]55]. For more details, see [115]Supplementary File S1.

   The poor yield (~30%) of the mapped reads with both genomes
   ([116]Tables S1 and S2) encouraged us to perform the transcriptome
   assembly of pure fish lines. First, trimmed reads were mapped against a
   custom rRNA database ([117]Supplementary File S2) using Bowtie2
   [[118]56] and also against the set of overrepresented sequences
   identified by FastQC using BLAST+ [[119]57]; in addition, all reads
   mapped to either rRNA or overrepresented sequences were removed to
   increase the protein-coding transcript yield. Filtered reads were
   checked using Rcorrector [[120]58] against erroneous k-mers, and all
   unfixable read-pairs were discarded.

   Transcriptome assembly was carried out by Trinity [[121]59], rnaSPAdes
   [[122]60] and MEGAHIT [[123]61] with multiple k-mer length values. All
   resulting transcripts were merged together using EvidentialGene’s
   tr2aacds [[124]62] on the basis of CDS-DNA local alignment identity
   classification. All transcriptome assemblies (including the merged one)
   were subjected to various quality control tools: TransRate [[125]63],
   rnaQUAST [[126]64] and BUSCO [[127]65]. Based on the quality assessment
   ([128]Tables S5–S8), we decided to proceed with the merged assembly.
   The annotation step involved TransDecoder [[129]66] and Trinotate
   [[130]67] as primary tools taking advantage of the following tools and
   databases: UniProtKB/Swiss-Prot database [[131]68], MEROPS database
   [[132]69], RefSeq database [[133]70], NCBI Nucleotide database
   [[134]71], BLAST+ [[135]57], Pfam database [[136]72], HMMER [[137]73]
   and SignalP [[138]74]. Finally, the abundance of the resulting
   transcripts was quantified using salmon [[139]75]. All compatible
   results and statistics were processed by MultiQC. For more details, see
   [140]Supplementary File S3.

   To be able to compare the expression data across lines, we would either
   need to cross-reference the two obtained transcriptomes or align all
   samples to the same reference. We tested the second option using the
   same approach we described earlier. Abramis brama uniquely aligned
   reads mapped to the R. rutilus transcriptome ended up only 4% less
   abundant (on average between samples) than when mapped to the A. brama
   transcriptome and vice versa ([141]Tables S3 and S4). As the R. rutilus
   transcriptome had more predicted genes assembled ([142]Tables S5 and
   S6), we decided to use it as our main reference and avoid the
   complicated cross-referencing process.

   Differential gene expression analysis between infected and non-infected
   individuals for each fish line was performed according to gene counts
   produced by featureCounts and further processed using the Bioconductor
   package DESeq2 [[143]76]. Genes were considered as differentially
   expressed only if having the adjusted p-value (FDR) smaller than 0.05
   and the absolute value of the log2 expression fold-change higher than
   1. The DESeq2 tool computes the normalized and shrunk log2 expression
   fold-change for each gene.

   The functional role of significantly differentially expressed genes
   (i.e., enrichment analysis) was performed according to the KEGG pathway
   database [[144]77] and the Gene ontology database [[145]78]. For each
   fish line, the differentially expressed genes were used as the target
   set, and all of the annotated genes within the transcriptome served as
   a background for the statistical testing separately for the up- and
   downregulated significantly differentially expressed genes. The Gene
   ontology (GO) annotation comes from the semi-automatic annotation made
   by Trinotate. From the total of 33,767 predicted genes, 12,590 were
   annotated with GO. KEGG annotation is derived from the KAAS—KEGG
   Automatic Annotation Server [[146]79] using GHOSTZ [[147]80] on peptide
   sequences of predicted genes against the set of all fish genes of the
   KEGG database. From the total of 33,767 predicted genes, 12,618 were
   annotated with KEGG orthology identifiers. Both functional analyses
   were done employing in-house scripts using the R packages
   clusterProfiler [[148]81] for the core of the analysis, data.table
   [[149]82] for postprocessing and making summary tables and pheatmap
   [[150]83] and ggplot2 [[151]84] for picturing the results
   ([152]Supplementary Material with Scripts S1 and S2).

   The expression change profile represents the difference in molecular
   response to P. homoion infection for a given gene. We calculated the
   Spearman correlation between the per-gene expression change profiles
   and the level of P. homoion infection. Biologically relevant genes with
   the highest correlation were selected for qPCR validation.

2.5. Gene Selection and Real-Time Quantitative PCR

   Three reference (housekeeping) genes, RPL8, 18S and α-tubulin (A-tub),
   were tested to normalize variation in the gene expression. To test
   their stability in our sample set, 67 representative samples (including
   at least 5 samples from each line) were subjected to analysis. The
   amplification was performed under the optimized conditions for qPCR
   mentioned below. The Reference Gene Selection Tool from Bio-Rad CFX
   Maestro software (Bio-Rad, Hercules, CA, USA), based on geNorm software
   principles [[153]85,[154]86] with an algorithm to normalize the Cq of
   each gene against the Cqs of all reference genes tested across the
   experimental samples (4), was used. On the basis of their stability,
   three tested references genes were recognized as follows: RPL8