Copy number increases of transposable elements and protein-coding genes in an invasive fish of hybrid origin.

Evolutionary dynamics of structural genetic variation in lineages of hybrid origin is not well explored, although structural mutations may increase in controlled hybrid crosses. We therefore tested whether structural variants accumulate in a fish of recent hybrid origin, invasive Cottus, relative to both parental species Cottus rhenanus and Cottus perifretum. Copy-number variation in exons of 10,979 genes was assessed using comparative genome hybridization arrays. Twelve genes showed significantly higher copy numbers in invasive Cottus compared to both parents. This coincided with increased expression for three genes related to vision, detoxification and muscle development, suggesting possible gene dosage effects. Copy number increases of putative transposons were assessed by comparative mapping of genomic DNA reads against a de novo assembly of 1,005 repetitive elements. In contrast to exons, copy number increases of repetitive elements were common (20.7%) in invasive Cottus, whereas decrease was very rare (0.01%). Among the increased repetitive elements, 53.8% occurred at higher numbers in C. perifretum compared to C. rhenanus, while only 1.4% were more abundant in C. rhenanus. This implies a biased mutational process that amplifies genetic material from one ancestor. To assess the frequency of de novo mutations through hybridization, we screened 64 laboratory-bred F2 offspring between the parental species for copy-number changes at five candidate loci. We found no evidence for new structural variants, indicating that they are too rare to be detected given our sampling scheme. Instead, they must have accumulated over more generations than we observed in a controlled cross.

Introduction

Hybridization can be an important driver of evolutionary change and may lead to homoploid hybrid speciation of animals and plants (Abbott et al., 2013; Mallet, 2007; Nolte & Tautz, 2010; Schumer, Rosenthal, & Andolfatto, 2014). It remains challenging to identify genetic factors that determine hybrid fitness, partly because of their complex genomic basis (Landry, Hartl, & Ranz, 2007). Reduced hybrid fitness can be explained through Bateson–Dobzhansky–Müller incompatibilities of alleles (Bateson, 1909; Dobzhansky, 1937; Müller, 1942) that can increase reproductive isolation and hamper the evolution of admixed lineages (Maheshwari & Barbash, 2011). On the other hand, hybridization can lead to new combinations of allelic variants and therefore be a source of evolutionary novelty and ecological niche divergence (Lexer, Lai, & Rieseberg, 2003). A common explanation is given by the process of transgressive segregation, which is caused by novel interactions of alleles that evolved in isolation and subsequently segregate in the hybrid offspring (Rieseberg, 1997; Rieseberg, Archer, & Wayne, 1999). This is plausible because hybridization introduces a wealth of genetic diversity that can affect evolutionary processes during the initial and highly dynamic phase of hybrid speciation.

In addition to effects caused by parental alleles, an increased rate of de novo structural mutations in hybrids appears possible, similar to the finding of increased point mutations in heterozygotes (“heterozygosity instability”; Amos, 2010; Xie et al., 2016). Whether this actually occurs and leaves detectable genetic signatures has received comparatively little attention in studies on hybrid speciation. However, it has already been found that hybrid crosses may exhibit de novo copy‐number variation that is frequent enough to be observed among siblings of the same litter of mice (Scavetta & Tautz, 2010). Although copy‐number variations (CNVs) may infer fitness costs associated with gene dosage alteration (Orozco et al., 2009; Papp, Pal, & Hurst, 2003), they may also contribute to the acquisition of new adaptive functions and promote species diversification (Kondrashov, 2012; Machado et al., 2014). Hybridization can also increase mutation rates through the activation of transposable elements (TEs), which have been reported to undergo “bursts of transposition” in hybrids of plants (McClintock, 1984; Renaut, Rowe, Ungerer, & Rieseberg, 2014; Ungerer, Strakosh, & Stimpson, 2009) and animals (García Guerreiro, 2015; Labrador & Fontdevila, 1994; O'Neill, O'Neill, & Graves, 1998; Renaut, Nolte, & Bernatchez, 2010). Transposition events can have numerous genomic consequences, including copy‐number changes of protein‐coding genes (Vogt et al., 2014; Xing et al., 2006), disruption of genes or regulatory motifs (Bonchev & Parisod, 2013) or the introduction of inversion break points that may negatively affect genomic stability (Casals, Cáceres, & Ruiz, 2003; Symer et al., 2002). Accordingly, increased transposition is often detrimental and therefore quickly silenced by the host genome (Castel & Martienssen, 2013). Conversely, transposons can also create adaptive genetic variation and be subject to positive natural selection (Casacuberta & González, 2013; Chuong, Elde, & Feschotte, 2016). Similar to copy‐number variation of genes, this has been repeatedly suggested to promote speciation and adaptive radiations (Belyayev, 2014; Fontdevila, 2005).

The expected number of gene and transposon copies in a lineage of hybrid origin should correspond with the average of the copy number in the parental species (e.g., Kawakami, Dhakal, Katterhenry, Heatherington, & Ungerer, 2011; Scavetta & Tautz, 2010). Deviations from this average would be conceivable if the genetic contributions from both parent species to the hybrid gene pool are not even. Copy numbers may, however, not exceed observed parental averages due to admixture alone. Instead, copy‐numbers outside the parental values could indicate the evolution of new copies and be a distinguishing feature of hybrid lineages. If rapidly evolving structural variants were also subject to positive selection, they could be among the first to contribute to adaptive evolution of hybrid lineages. Although the underlying mutational mechanisms are known from controlled hybrid crosses (e.g., García Guerreiro, 2014; McClintock, 1984; Scavetta & Tautz, 2010), few empirical studies have documented the evolutionary dynamics of copy‐number variation and transposition in natural admixed lineages. For example, studies on homoploid sunflower hybrid species revealed a massive accumulation of Ty3/gypsy‐like TEs during the past 0.5–1 million years (Renaut et al., 2014; Ungerer et al., 2009). On the other hand, transposition bursts have not been observed in F1 hybrid sunflowers or contemporary admixture zones between their parental species (Kawakami et al., 2011; Renaut et al., 2014). Consequently, the relative importance of short‐term compared to long‐term evolutionary processes that drive the proliferation of TEs remains unclear when hybrid lineages are ancient.

Here, we analysed whether CNVs have increased in a small freshwater fish of hybrid origin that has emerged postglacially, possibly within the last 200 years (Nolte, Freyhof, Stemshorn, & Tautz, 2005). The so‐called invasive Cottus provides a unique natural setting to explore evolutionary mechanisms that may initiate the evolution of a homoploid hybrid species (Nolte & Tautz, 2010). Invasive Cottus is a lineage of hybrid origin that has colonized lower reaches of rivers and streams of intermediate size, whereas the parental species Cottus perifretum and Cottus rhenanus live in parapatry, apparently confined to the headwaters of smaller tributaries (Nolte et al., 2005; Stemshorn, Reed, Nolte, & Tautz, 2011). A likely explanation for the success in lower reaches of rivers is that invasive Cottus has evolved new beneficial traits through hybridization that are driving divergence between invasive Cottus and its parental species. Besides barriers to gene flow in hybrid zones (Nolte, Freyhof, & Tautz, 2006; Nolte, Gompert, & Buerkle, 2009), this manifests in patterns of gene expression that are unique to invasive Cottus (Czypionka, Cheng, Pozhitkov, & Nolte, 2012). In this study, we tested whether rapid copy‐number evolution of protein‐coding genes and TEs occurs in natural populations of invasive Cottus. We expected to find overall more gene duplications than deletions (e.g., Katju & Bergthorsson, 2013), as well as an increased variance of copy numbers compared to the parental species (Scavetta & Tautz, 2010). To further explore the possibility that hybridization causes rapid evolution of de novo structural variants, we tested laboratory‐raised F2 (selfed F1) C. rhenanus × C. perifretum crosses for copy‐number changes of candidate genes and transposons that have increased in copy number in natural populations of invasive Cottus. Finally, we searched for associations between changes in gene copy number and gene expression in invasive Cottus to explore possible contributions of CNVs to phenotypic change.

Materials and Methods

Sample collection and DNA extraction

Samples were collected in Northern Germany (watersheds of rivers Sieg, Wied, Emscher, Mosel, Lippe, Ruhr, Wupper, Wisper), Belgium (River Scheldt watershed) and Great Britain (watersheds of rivers Wensum, Great Ouse) between 2001 and 2015 (Table 1). Gender was determined by carefully inspecting gonad development. Genomic DNA was extracted from blood samples, using the 5 prime ArchivePure DNA Tissue kit (5 PRIME, Hamburg, Germany) or the Qiagen Blood and Tissue kit. Extracted DNA was quality checked for the presence of high molecular weight genomic DNA by agarose gel electrophoresis.

Table 1

Individual samples used for comparative genomic hybridization array (aCGH) microarrays and genome mapping. Indicated are number, gender and geographic coordinates for each population

Population	aCGH	Genome mapping	Lat./Long.
Cottus rhenanus
Bröl (Sieg)	1fa	1	50°47′43″N 7°20′22″E
Mengbach (Sieg)	2f	1f	50°46′57″N 7°25′13″E
Naafbach (Sieg)	1ma	1	50°52′00″N 7°16′00″E
Krabach (Sieg)	1m, 1f	1m	50°45′45″N 7°24′47″E
Wahnbach (Sieg)	1m	1m	50°50′50″N 7°19′15″E
Ottersbach (Sieg)		1m	50°46′32″N 7°28′56″E
Fockenbach (Wied)	1m	1	50°32′14″N 7°25′46″E
Boye (Emscher)		1m	51°33′33″N 6°55′42″E
Kyll (Mosel)		1	50°14′00″N 6°39′00″E
Soestbach (Lippe)		1m	51°34′45″N 8°05′17″E
Wisper (Rhein)		1	50°20′58″N 7°50′60″E
Wupper (Rhein)		1	51°15′22″N 7°16′19″E
Wanne (Ruhr)		1m	51°25′29″N 8°03′05″E
Cottus perifretum
Larse Beek (Schelde)	2fa	1	51°17′10″N 4°30′11″E
Witte Nete (Schelde)	3ma, 1fa	1	51°14′22″N 5°4′11″E
Zwanebeek (Schelde)	2	1	51°14′44″N 4°34′58″E
Trouille (Schelde)		1	50°23′55″N 4°0′52″E
Molenbeek (Schelde)		1	50°48′14″N 3°40′42″E
Twin (Great Ouse)		1	51°58′10″N 0°57′59″E
Wensum (Yare)		1	52°42′59″N 1°1′59″E
Cottus “invasive”
Sieg (Rhein)	4m, 4f	5m, 2f, 3	50°47′55″N 7°10′45″E

Stream name with higher order river given in parentheses.

m = male, f = female.

Laboratory‐bred.

Array‐CGH design and calibration

We developed a custom 60‐mer oligonucleotide comparative genomic hybridization array (aCGH) to screen for copy number differences across 10,979 genes that could be annotated to the threespine stickleback genome (Gasterosteus aculeatus; BROADS1.56). Briefly, an aCGH experiment compares genomic DNA of test and reference samples that are digested (optional), colour‐labelled (Cy3; Cy5) and hybridized together on an array that is scanned for fluorescence intensity values. The presence of gene duplications or deletions is inferred from log2‐intensity ratios between test and reference signals.

Probe design for the microarrays was largely based on Cottus expressed sequences (ESTs) assembled with the trinity software (Grabherr et al., 2011). For the transcriptome assembly, total RNA was extracted from different tissues (gills/heart, liver, fins/skin, digestive tract, head without gills, white muscle) of a C. rhenanus male from Fockenbach using a TRIzol protocol (Invitrogen, Carlsbad, CA, USA). Libraries were prepared with the Illumina Truseq Kit and sequenced from two sides (100‐bp read length) on an Illumina GAIIx sequencer. We only targeted genes that could be annotated to the threespine stickleback genome, aiming for a coverage of one probe per exon (“exon‐array CGH”; Dhami et al., 2005). Additionally, we included probes for genes that were identified to be either upregulated or downregulated in invasive Cottus in a previous gene expression microarray study (Czypionka et al., 2012). We first conducted a blat search (blat version 35x1; Kent, 2002; e‐value < 10−3) of Cottus ESTs against the stickleback genome and subsequently used best matching sequences larger than 60 bp for a second blat search against stickleback coding sequences (retrieved from the biomart database version 0.7; www.ensembl.org/biomart). The sense direction of all probe sequences was inferred from the directionality given in the blat results. Sequences identified as negative strands were reverse‐transcribed, and sequences with contrasting directionality information between the ensemble database and the blat results were removed.

To optimize probe selection and to validate the performance of probes, we initially performed a microarray calibration experiment (Czypionka et al., 2012; Pozhitkov, Noble, Bryk, & Tautz, 2014). We designed three probes per exon, which were selected according to an optimal melting temperature of 80°C. Probe sequences were submitted to SureDesign (www.agilent.com/genomics/suredesign) for the design of a SurePrint G3 custom CGH microarray (8x60k) with two replicates per probe and the addition of a 15‐bp linker recommended by Agilent. A dilution series of eight starting DNA concentrations ranging between 0.0625× and 4× of a recommended starting DNA amount of 400 ng (corresponding to 25, 50, 100, 200, 400, 800, 1,200, and 1,600 ng) was prepared for a single C. rhenanus and a single C. perifretum individual. Both dilution series (C. rhenanus, Cy3 labelled; C. perifretum, Cy5 labelled) were cohybridized on the calibration array. Linear regressions on average (of two probe replicates) signal intensities across the dilution gradient were fitted for red (Cy5) and green (Cy3) channels separately. Next, we excluded probes with nonlinear response curves (R 2 < .95), a signal‐to‐error ratio below 10, as well as unresponsive probes. From the remaining probe set, we retained one probe per gene by filtering for a signal intensity increase (regression slope) being closest to one. The final probe set contained 12,105 probes covering 10,979 stickleback genes, as well as 446 candidate probes from Czypionka et al. (2012) (see Table S3). The final probe set was sent to SureDesign for the construction of an experimental 8x60k array with five replicates per probe, and a 15‐bp linker added to each probe as recommended by Agilent (SureDesign‐ID: 070467).

The experimental array was used for a common reference aCGH experiment. Three array slides were used to cohybridize 24 test samples (eight individuals from each of C. rhenanus, C. perifretum, and invasive Cottus; all Cy5 labelled) with a common reference sample (a female individual of C. perifretum from Larse Beek; Cy3 labelled) (Table 1). Samples were randomly assigned to the 24 arrays, and for each sample, genomic DNA was quantified using a fluorescence NanoDrop (dsDNA BR Assay Kit, Thermo Fisher Scientific, Waltham, MA, USA) to achieve equal amounts (400 ng) of starting DNA. For all aCGH experiments, genomic DNA was labelled following the “direct method” as outlined in the Agilent protocol “Agilent Oligonucleotide Array‐based CGH for genomic DNA analysis” (version 7.3, March 2014; Agilent Technologies). The same protocol was used for subsequent hybridization and washing of the arrays. Microarrays were covered with Agilent ozone‐barrier slide covers and scanned using an Agilent Microarray Scanner (type C). Signal intensity values were extracted using the Agilent feature extraction software version 10.7.3.1.

Quality filtering, normalization

Quality control metrics provided in the feature extraction output files were carefully inspected for the spatial distribution of outliers on the arrays. A pipetting error before hybridization created a noticeable bubble in the centre of nine experimental arrays, resulting in decreased signal intensity values and clustering of outliers in these regions. We therefore decided to remove all data affected by such bubble artefacts (~12% of 1,452,600 probe replicates present on 24 arrays). For this purpose, we first retrieved coordinates for individual probe replicates (“features”) from the feature extraction output files and plotted the spatial distribution of log‐signal intensities for each array (using the “levelplot” function in the r package latticeextra). We used the r package “splancs” (version 2.01‐36; Bivand & Gebhardt, 2000) to extract the coordinates of features located within bubble artefacts by drawing polygons around the bubbles, which were clearly visible through reduced log‐signal intensities. Coordinates of all features located within bubble areas were then used to remove affected data from the complete data set. Next, we excluded outliers among probe replicates as determined by signal intensities exceeding the 1.5× interquartile range and took the average of the remaining probe replicates for further downstream analyses. The final data set covered 12,102 probes represented on average on 22 out of 24 arrays.

To account for probe‐specific binding coefficients in the observed signal intensities, we applied a calibration on average signal intensities as recommended by Pozhitkov et al. (2014). For this purpose, we used the probe‐specific linear regression parameters (slope, intercept) from the calibration array, to estimate a calibrated starting DNA value (DNAcal) for each channel (Cy3, Cy5), calculated as DNAcal = (observed signal − intercept)/slope. Although this procedure reduces artefacts introduced through probe‐specific differences in binding behaviour, it does not remove intensity effects, which are visible in a “MA”‐plot as increasing log2‐(Cy5/Cy3)‐ratios (“M”) with increasing average (log‐) signal intensities (“A”). To remove such intensity effects and to normalize log2‐ratios to an average of zero, we performed a LOESS normalization (span = 0.3, iterations = 4) without background correction (Zahurak et al., 2007), as implemented in the r package limma (version 3.22.4, Ritchie et al., 2015). This method fits a smoothing curve based on locally weighted least‐squares regression and uses the residuals of the fitted curve as normalized log2‐ratios.

Array‐CGH data analyses and comparison to gene expression

Genes showing increased or decreased copy‐numbers in invasive Cottus compared to the (combined) parental species were identified using two‐sample Wilcoxon rank‐sum tests on normalized signal intensity values in r (FDR corrected for multiple testing; Benjamini & Hochberg, 1995). To test whether invasive Cottus show a higher tendency towards increased compared to decreased copy numbers, we tested for skewness of average log2‐ratios of invasive Cottus relative to parental Cottus log2‐ratios using a two‐sided D'Agostini test of skewness implemented in the r package “moments.” To test whether invasive Cottus showed an overall increase in variability of copy numbers, we conducted t tests (FDR corrected for multiple testing) to assess how many probes showed significantly higher variance of log2‐ratios in invasive compared to parental Cottus, and how often each of the parental species showed significantly higher variance compared to the other species.

Possible gene dosage effects of exon copy‐number changes in invasive Cottus were evaluated by comparing CNV candidates with lists of gene expression candidates showing significant higher expression in invasive Cottus compared to the parental species. These gene expression candidates were characterized in a previous microarray study using whole‐body tissues of juvenile Cottus individuals (Czypionka et al., 2012).

Genome sequencing, Repetitive element assembly and CNV analyses

We investigated the possibility of transposon accumulation in invasive Cottus by comparing the relative number of genomic reads that could be mapped to a de novo assembly of repetitive elements between invasive and parental Cottus species. For this purpose, we used whole‐genome Illumina paired‐end sequence reads of 13 C. rhenanus, seven C. perifretum and 10 invasive Cottus genomes (Table 1). Genomic DNA was fragmented using sonication (Covaris, Woburn, MA, USA), and fragments were end‐repaired and adapter ligated. The libraries were sequenced using an Illumina GAIIx sequencer with 2× 150‐bp paired‐end reads. Raw paired‐end Illumina GAII reads were quality trimmed using the software popoolation (Kofler et al., 2011). To facilitate the assembly process, sequence reads were first error‐corrected using the software lighter, with a k‐mer length of 31, a genome length of 1 Gbp, and a subsampling fraction of 0.1 (Song, Florea, & Langmead, 2014). For each individual, we created assemblies of putative TEs using a de Bruijn graph‐based approach implemented in the software tedna v. 1.2.2 (Zytnicki, Akhunov, & Quesneville, 2014). As parameter settings, we used a maximum graph size of 1,000, kmer 30, minimum size of 500 bp, and an insert size of 300 bp. To create an improved assembly with longer contigs, we assembled the resulting tedna output sequences in the software seqman pro version 2.0.0 (DNASTAR, Madison, WI, USA) using default settings. The consensus sequences were exported as fasta files and retained as a Cottus reference library of repetitive elements. To classify putative transposons, we used the software pasteclassifier (v. 1.0; Hoede et al., 2014). pasteclassifier searches for characteristic structural and functional features of TEs, including long terminal repeats (LTR), terminal inverted repeats (TIR), reverse transcriptases (RT), transposases, polyA tails and open‐reading frames (ORFs), and tests for similarities to known TEs provided by the repbase database (v. 20.04; Jurka et al., 2005), as well as to hidden Markov model‐based profiles of complete protein domains (pfam26.0_gypsydb; Finn et al., 2016). By combining all this information, pasteclassifier then classifies putative transposons to class and order level according to the Wicker hierarchical TE classification system (Wicker et al., 2007). Remaining unclassified sequences were further analysed using blastx searches in blast2go (Conesa et al., 2005).

To compare CNV of putative TEs among Cottus individuals, we mapped whole‐genome (Illumina) sequence reads (trimmed to 90‐bp length to ensure comparability among genomes) against the Cottus reference library of repetitive elements using the software nextgenmap with a minimum identity constraint of 95% (ngm version 0.4.12, Sedlazeck, Rescheneder, & von Haeseler, 2013). For each Cottus individual, we counted the number of mapped forward reads, divided it by the total number of reads that were used for mapping and multiplied it by one million. One‐sample Wilcoxon rank‐sum tests (FDR corrected for multiple testing) were then used to compare counts per million mapped reads between invasive and parental Cottus species.

Validation of candidate CNVs and screen for de novo CNVs in F2 crosses

We used digital droplet PCR (ddPCR; upgraded QX100 Droplet Digital PCR System; Bio‐Rad Laboratories, CA, USA) with EvaGreen dye as fluorescent marker to (1) validate a copy‐number increase of three genes and three repetitive elements for the same 24 individuals used for array‐CGH experiments and to (2) screen laboratory F crosses for de novo structural variants of ddPCR‐validated candidate CNVs. For the de novo structural variant screen in F crosses, we used 32 individuals from each of two C. rhenanus × C. perifretum families that were obtained by selfing F1 hybrids between individuals representing the ancestral species (Populations crossed: Bröl × Witte Nete and Naaf × Larse Beek, respectively).

To estimate copy numbers, we used multiplex reactions of the target sequence together with a single‐copy reference gene (McDermott et al., 2013; Miotke, Lau, Rumma, & Ji, 2014). As a reference single‐copy gene, we chose rpl13a, which has been used as a single‐copy housekeeping gene in threespine stickleback (Hibbeler, Scharsack, & Becker, 2008), and showed no signs of CNV in our own CGH‐array experiments. All primers were designed using primer3 (http://primer3.ut.ee) and together with the targeted regions visualized in igv Viewer (version 2.3.32, Thorvaldsdóttir, Robinson, & Mesirov, 2013) to ensure coverage by all individuals and to avoid common SNPs. Initial gradient PCRs for annealing temperature and different DNA concentrations were performed to optimize PCRs. For all PCRs, 20 μl of the mastermix was mixed with 70 μl Droplet Generation Oil (Bio‐Rad Laboratories) and partitioned into droplets using a QX100 Droplet Generator (Bio‐Rad Laboratories). For primer sequences as well as ddPCR cycling conditions, see Table S1. For all PCR products, fluorescence of individual droplets was measured using an upgraded QX100 Droplet Digital reader (Bio‐Rad Laboratories) and analysed with the quantasoft droplet reader software (Bio‐Rad Laboratories). To test for significant differences of copy numbers among species, we applied one‐way ANOVAs followed by Tukey HSD post hoc tests.

Results

Exon‐Array CGH copy‐number variation

We compared CNV in eight invasive Cottus individuals with eight individuals of each of the parental species, C. rhenanus and C. perifretum to detect copy‐number changes of genes in invasive Cottus. Among 10,979 genes (exons) and additional 446 probes derived from a previous gene expression study (Czypionka et al., 2012), we identified 12 genes with increased as well as 13 genes with decreased copy numbers in invasive Cottus compared to both parental species (two‐tailed Wilcoxon rank‐sum test; FDR corrected; p < .05) (Figure 1; Table S2). Between the parental species C. perifretum and C. rhenanus, we found 1,920 genes (~15.9%) to differ in copy number (two‐tailed Wilcoxon test; FDR corrected; p < .05). Copy numbers did not differ significantly between males and females (results not shown).

Figure 1

Candidate copy‐number variations genes from aCGH analyses with significantly increased (a) and decreased (b) log2‐ratios in invasive Cottus (dark grey = invasive Cottus, light grey = C. perifretum, white = C. rhenanus)

In accordance with the general expectation to find more gene duplications than deletions, copy‐number variable genes in invasive Cottus showed a stronger signal of copy number increase than decrease (D'Agostini test: skew = 0.8918; z = 34.6861; p < .001). This trend of invasive Cottus to show more increases of copy numbers was also visible using arbitrary log2‐thresholds: Increased copy numbers were inferred for 409 genes that exceeded a log2‐ratio threshold of 0.5 (putative duplications), whereas decreased copy numbers were inferred for 314 genes falling below a log2‐ratio threshold of −0.5 (putative deletions) in at least one individual. Among those genes, 125 included individuals with signs of both deletions and duplications, leaving a total of 598 copy‐number variable genes (4.9%). When applying a less stringent log2‐ratio threshold of (±) 0.4, a total of 1346 (11.1%) of all genes were estimated to be copy number variable in at least one individual, which included 916 putative duplications and 733 putative deletions.

In contrast to our expectation, we found no increased variance of copy numbers in invasive Cottus compared to the parental species. Compared to C. rhenanus, the number of genes showing an increased variance was 825 in invasive Cottus and 947 in C. perifretum. Cottus rhenanus showed similar numbers of genes with increased variance compared to both other species: 769 genes were more variable than in C. perifretum and 685 genes more variable than in invasive Cottus. The variances between invasive Cottus and C. perifretum differed conspicuously from the above numbers in that differences between genes were comparably rare: invasive Cottus showed a higher variance than C. perifretum in 33 genes, whereas C. perifretum showed a higher variance than invasive Cottus in 36 genes.

Comparison of CNV candidates to gene expression

To explore possible gene dosage effects of increased exon copy numbers, we compared CNV candidate genes to lists of genes that were overexpressed in invasive Cottus. Possible gene dosage candidates included four probes (GTF2IRD2, CYP27C1 and two probes for CYP2K1) that were previously found to be overexpressed in whole‐body tissues of juvenile individuals of invasive Cottus (Czypionka et al., 2012).

De novo assembly and copy‐number variation of repetitive elements

To assess a possible accumulation of TEs in invasive Cottus, we used a genome mapping approach to compare the relative number of reads mapped to a de novo assembly of repetitive elements. The de novo assemblies in tedna resulted in a total of 7,548 sequences, which were assembled in seqman pro into 1,005 contigs (hereafter referred to as “repetitive elements”; size range: 500‐50k bp; n50: 2,978 bp; Table S4). Using pasteclassifier, we classified 452 repetitive elements (45%) according to Wicker's classification system of TEs. These included 312 putative retrotransposons (Class I elements), 119 putative DNA transposons (Class II elements) and 21 other elements (potential host genes). Class I transposons included 110 LTR, 7 DIRS, 8 PLE, 6 LARD, 20 TRIM, 137 LINE, 11 SINE and 13 undetermined elements. Class II transposons included 87 TIR, 3 MITE, 12 Crypton, 5 Helitron, 1 Maverick and 11 undetermined elements. Potential host genes included for example a large fraction of the mitochondrial genome, various proteins, zinc finger domains, RNA binding and recognition motifs, and a herpesvirus tegument protein. Among the unclassified 553 repetitive elements, a blastx search implemented in blast2go matched 180 elements to various host gene proteins, one herpesvirus‐homolog, one transcription factor IIIA‐like element (GTF3A) as well as 12 retrovirus‐like polyproteins (pol; gag‐pol).

To assess whether invasive Cottus individuals (n = 10) showed higher copy numbers compared to both parental species (n = 20), we applied Wilcoxon rank‐sum tests to compare individual count data for each of 1,005 repetitive elements (one‐tailed; FDR corrected). Among all 1,005 repetitive elements, we found significant copy‐number increases in 208 (20.7%) repetitive elements with an average copy‐number increase of 30.5% compared to the parental mean (Figure 2; marked by red sidebar). In contrast, significant decreases in copy number were rare, affecting only 12 (0.01%) repetitive elements (nine unclassified elements, two putative retrotransposons, one potential host gene). Overall, increased repetitive elements showed a strong bias towards those being present in higher copy numbers in C. perifretum compared to C. rhenanus (Figure 2). This was also evident when comparing copy numbers between C. perifretum and C. rhenanus for the 208 candidate elements, revealing 112 elements (53.8%) with significantly higher copy numbers in C. perifretum and three elements (1.4%) with significantly higher copy numbers in C. rhenanus (one‐tailed Wilcoxon rank‐sum tests; FDR corrected; p < .05).

Figure 2

Heatmap showing mapping results of genomic reads from 30 Cottus genomes against 1,005 repetitive elements. Colour‐coding represents the percentage of sum of reads (number of mapped reads per Mio reads) for each repetitive element. Coloured sidebars on the left indicate 208 repetitive elements with significantly higher numbers of mapped reads in invasive Cottus compared to the parental species (red), 12 elements with significantly lower numbers in invasive Cottus than in the parental species (orange), as well as repetitive elements that are significantly higher in C. rhenanus (light grey), and significantly lower in C. rhenanus (dark grey) compared to combined invasive Cottus and C. perifretum

Out of the 208 candidate elements, 86 (41.3%) could be classified using pasteclassifier. These included 76 putative retrotransposons (Class I elements: 18 LTR, 2 DIRS, 4 PLE, 1 LARD, 7 TRIM, 38 LINE, 3 SINE), 19 putative DNA transposons (Class II elements: 12 TIR, 1 MITE, 2 Crypton elements) and 1 potential host gene (EFHD2). From the remaining 122 unclassified candidate elements, 27 showed blastx matches to various proteins such as GTF3A (induces transcription of 5S rRNA), interferon‐induced 44 protein (a putative response to hepatitis virus) and G2E3 (a ubiquitin‐protein ligase).

Among the 86 classified candidate transposons, 15 were marked as “complete” by the pasteclassifier software, which included mostly non‐LTR retrotransposons (11 LINE, 1 SINE, 2 PLE, 1 TIR elements). No significant overrepresentation of either Class I or Class II elements could be detected among candidate transposons that were increased in invasive Cottus (Fisher's exact test; two‐sided; p = .264). When subdividing Class I elements further into LTR‐ and non‐LTR retrotransposons, elements classified to order showed a nonsignificant (Fisher's exact test; two‐sided; p = .053) bias towards non‐LTR retrotransposons.

Overall, the proportion of increased elements in invasive Cottus (20.7%) was similar between classified transposons and unclassified elements. Of 452 repetitive elements classified to transposon‐order, 86 (19%) were increased in invasive Cottus (retrotransposons: 24.4%; DNA transposons: 16%). Among 553 unclassified repetitive elements, 122 (22.1%) were increased in invasive Cottus, and out of 180 unclassified repetitive elements with some blastx match; 27 (15%) were increased in invasive Cottus.

Validation of CNV candidates and screen for de novo CNVs in F2 crosses

We used ddPCR to validate increased copy numbers in invasive Cottus for three selected CNV candidate genes (MIEF1, NUGGC, CYP27C1) and three repetitive elements (one LINE element, two uncharacterized elements). These CNV candidates showed some of the strongest signals of copy‐number increases in invasive Cottus and were therefore chosen as a proof of principle to confirm our alternative approaches to detecting copy‐number variants. In the parental species, ddPCR revealed CYP27C1 as a single‐copy gene, whereas MIEF1 occurred mostly in two copies and NUGGC was present as a more variable multicopy gene with most individuals showing five to seven copies (Figure 3).

Figure 3

Digital droplet PCR (ddPCR) validations of copy‐number variations candidates for three genes (a–c) and two repetitive elements (d, e). Box plots show copy‐number estimates for invasive Cottus (INV), C. perifretum (PER), C. rhenanus (RHEN) as well as two F2 offspring families from laboratory crosses between C. rhenanus and C. perifretum (BxWN = Bröl × Witte Nete; NxLB = Naaf × Larse Beek)

Because the CYP27C1 gene expression candidate probe from Czypionka et al. (2012) was located outside exonic regions, we additionally tested a newly designed probe located within (putative) exon 8 to test for a duplication of the protein‐coding gene. Increased copy numbers could be confirmed (ANOVA; followed by Tukey post hoc tests; p < .05) for the three genes (including both exonic and nonexonic probes of CYP27C1), the LINE element and one uncharacterized repetitive element (Figure 3). One uncharacterized element (“Contig_926”) showed similar high copy numbers in invasive Cottus and C. perifretum. Among two of the increased genes (MIEF1, NUGGC) and the two increased repetitive elements, single individuals of either of the parental species showed similarly high copy numbers as observed in invasive Cottus (Figure 3).

Candidates with strong signals in our analyses were considered useful to test for de novo mutations in controlled crosses. To explore a possible role of rapidly evolving de novo structural variants, we screened 32 individuals from each of two F2 families (C. rhenanus × C. perifretum) for CNV of validated candidate genes (MIEF1, NUGGC, CYP27C1) and repetitive elements (one LINE element, one uncharacterized element) (Figure 3). None of the examined CNV candidates showed signs of de novo copy‐number increases compared to the parental species (Figure 3).

Discussion

Copy‐number variation has long been recognized to play a role in the evolution of phenotypic novelty. However, the possible contribution of structural variations to the evolution of admixed lineages remains largely unexplored in natural systems. Here, we report copy‐number increases of candidate genes (exons) as well as a proliferation of repetitive elements in a homoploid hybrid lineage of Cottus. The origin of invasive Cottus is most likely postglacial and the biogeographic settings suggest an age of <200 generations (Nolte et al., 2005). Hence, copy number increases that are specific to invasive Cottus could have evolved within a relatively short time frame following admixture. Herein, we found no evidence for de novo structural mutations in three genes and two repetitive elements in two families of laboratory F2 crosses between the parental species. However, this limited set of 64 F2 offspring only represents an upper bound for the frequency of such mutations in early generation hybrid crosses. It is still possible that de novo mutations occurred in natural hybrid populations at a frequency that could not be detected here. On the other hand, our ability to distinguish whether CNVs in admixed lineages are recruited from standing genetic variation or from de novo mutations is limited because of uncertainty about the variation present in natural populations. We discuss evolutionary processes that may be relevant to explain copy‐number changes in invasive Cottus and how they could be involved in adaptive evolution.

Exon copy‐number variation in admixed lineages

Although copy‐number variation between the parental species occurs frequently, we found only a very low proportion of copy number variable exons that may have increased in invasive Cottus. This is not surprising given the presumably recent origin of invasive Cottus and considering that gene duplication rates are typically slow (Katju & Bergthorsson, 2013; Lynch & Conery, 2000). It is also consistent with the general assumption that new gene duplications are quickly removed from the genome by purifying selection (Katju & Bergthorsson, 2013; Lynch & Conery, 2000; Ohno, 1970). However, the occurrence of gene duplicates that reach high frequencies in invasive Cottus is conspicuous. Genomic admixture in hybrids can increase structural mutation rates through retrotransposition (Michalak, 2009) or nonallelic homologous recombination events, which can be promoted by decreased genomic stability particularly in regions with a high density of tandemly arrayed repeats (Janoušek, Karn, & Laukaitis, 2013). Accordingly, copy number increases of some genes that are present in multiple copies in the parental species could be explained through nonallelic homologous recombination between tandemly arrayed gene copies. This could be true for NUGGC, which occurs in variable multicopy numbers in the parental Cottus species. A tandem organization of NUGGC copies appears likely because this can also be observed in the genome of the closely related threespine stickleback (UCSC Genome Browser). A larger segmental duplication could have played a role in the duplication of MIEF1. For this gene, conservation of synteny between Cottus and stickleback (Cheng, Czypionka, & Nolte, 2013) infers si:ch211‐256m1.8 as the most likely neighbouring gene, which shows a joint increase of copy number in invasive Cottus. Increased copy numbers could also originate from standing genetic variation already present in populations or individuals of the parental species not sampled here. This seems plausible for genes such as MIEF1 or NUGGC, for which ddPCR identified single individuals of the parental species with increased copy numbers. Hence, our study did not yield conclusive evidence that genomic admixture caused de novo copy number increases of coding genes in invasive Cottus.

While we cannot rule out genetic drift, the occurrence of high‐frequency gene duplicates that are close to fixation in invasive Cottus but not the parental species indicates a possible role of selection. A recent effort to disentangle the relative role of drift and selection suggested that about half of the high‐frequency gene duplicates in Drosophila melanogaster were fixed by selection and the other half by drift (Cardoso‐Moreira et al., 2016). Assuming lower effective population sizes in Cottus, this could mean a higher proportion of gene duplicates being fixed by drift. However, the time to fixation should be 4Ne under neutral expectations (e.g., Innan & Kondrashov, 2010), which in the case of 200 generations assumed for invasive Cottus would translate into a Ne of 50 that seems unrealistically low. Drift is therefore unlikely to explain all gene duplicates, and instead, some of the gene duplicates may have risen to high frequencies through selection (e.g., Feulner et al., 2013; Kondrashov, 2012). Selection could favour new alternative splice variants following duplications of exons (Abascal, Tress, & Valencia, 2015) or act on gene dosage effects that are thought to underlie most of the reported adaptive gene duplicates (Kondrashov, 2012). Increased gene dosage appears to be a possible consequence of duplications in a GTF2IRD2‐like gene and two members of the cytochrome P450 gene family (CYP27C1, CYP2K1). These genes are among the strongest candidates for increased copy numbers in our aCGH analyses and were overexpressed in juvenile stages of invasive Cottus in a previous gene expression study (Czypionka et al., 2012). Functions associated with cytochrome P450 genes are highly diverse (Danielson, 2002) and have been related for example to insecticide resistance (Wondji et al., 2009), adaptation to desert environments (Jirimutu et al., 2012) or immune responses to bacterial and virus infections (Iizuka et al., 2004). In teleost fishes, CYP2K1 appears to be mostly expressed in liver and kidney and involved in detoxification of environmental pollutants (e.g., Buhler, Zhao, Yang, Miranda, & Buhler, 1997; Uno, Ishizuka, & Itakura, 2012). Enright et al. (2015) reported that CYP27C1 may red‐shift the spectral sensitivity in zebrafish by converting vitamin A1 into A2, which translates into an enhanced behavioural response to long‐wavelength light. If CYP27C1 has a similar function in Cottus, increased gene dosage could be beneficial for enhanced vision in more turbid environments (Fain, 2015) such as the larger streams inhabited by invasive Cottus. Finally, GTF2IRD2 might affect the development of skeletal muscle during the growth period of juvenile stages, with increased gene dosage leading to a higher proportion of slow‐type muscle fibres (Palmer et al., 2012). Changed proportions of muscle fibre types could be involved in the acclimation to different temperature regimes by invasive Cottus, resembling what has been reported for other teleost fishes (e.g., Woytanowski & Coughlin, 2013). Altogether, we found some gene duplicates with potentially adaptive gene dosage effects in an evolutionary young homoploid hybrid species, but their low number and the possibility that they were recruited from standing genetic variation of the parental species does not suggest a major role of increased de novo evolution of structural variants.

Increases of repetitive elements

We detected a significant increase of copy numbers in 20.7% of repetitive elements including putative transposons in invasive Cottus. The frequency of these copy‐number increases exceeded what we observed with protein‐coding genes by far. Their proliferation could have been be induced by a “genomic shock” following exposure to environmental stress or hybridization (Belyayev, 2014; Fontdevila, 2005; McClintock, 1984; Senerchia, Felber, & Parisod, 2015). A weak bias towards non‐LTR retrotransposons among increased TEs is in line with such a burst of transposition, but this excess is not significant and other repetitive elements have increased almost equally. For example, we find increased copy numbers of putative DNA transposons that use a cut‐and‐paste mechanism and should only rarely increase in copy number. Similarly, a range of uncharacterized repetitive elements that are not likely to multiply via transposition also increased. Hence, the accumulation of repetitive elements in invasive Cottus cannot be fully explained through transposition. Resembling our results, transposition bursts are not evident in F1 hybrid sunflowers or contemporary hybrid zones between their parental species (Kawakami et al., 2011; Renaut et al., 2014). This has led to argue that a combination of hybridization and environmental stress during colonization of extreme habitats (deserts, salt marshes) could have initiated bursts of transposition in hybrid sunflower species (Kawakami et al., 2011), although subsequent experimental tests failed to induce TE proliferation under various stress conditions (Ungerer & Kawakami, 2013). The absence of transposition bursts has also been reported from retrotransposons in polyploid plant hybrids, which instead emphasized the importance of TE silencing through reorganization of methylation states and deletions of TE loci (e.g., Senerchia et al., 2015, 2016). Accordingly, an accumulation of TEs through transposition activities may be rare and is presumably limited to only a small fraction of TEs that escape silencing mechanisms without provoking strong deleterious effects (Parisod et al., 2010).

As an underexplored alternative, increased copy numbers of repetitive elements could be explained through multiple segmental duplications in centromeric or pericentromeric regions, which typically harbour large arrays of satellite repeats and TEs (Horvath et al., 2005; Ma & Jackson, 2006). Centromeres have been suggested as a primary target for genomic destabilization following interspecific hybridization in marsupial hybrids, which was shown to cause amplification of satellite and TEs (Metcalfe et al., 2007). Likewise, a massive accumulation of retro TEs has occurred in pericentromeric regions of homoploid hybrid species of sunflowers (Staton, Ungerer, & Moore, 2009). Interestingly, most of the repetitive elements increased in invasive Cottus have higher copy numbers in C. perifretum compared to C. rhenanus (Figure 3). Conversely, repetitive elements that occur in higher copy numbers in C. rhenanus compared to C. perifretum have rarely expanded in the invasive Cottus genome. This asymmetric pattern of repetitive element amplification correlates with an overrepresentation of C. perifretum ancestry in the genome of invasive Cottus, which contributes ~60% to the hybrid gene pool (Stemshorn et al., 2011). Thus, a nonrandom mechanism of repetitive element amplification with respect to ancestry appears possible. Centromere drive could quickly create such an asymmetry though a transmission advantage of longer centromeres during female meiosis, mediated by a higher number of microtubule binding sites during the segregation of chromosomes (Henikoff, Ahmad, & Malik, 2001). Such a transmission bias could produce an asymmetric pattern of repetitive element accumulation if amplified elements in invasive Cottus are indeed located in centromeric regions. If true, enlarged centromeres from C. perifretum but not C. rhenanus could sweep through the hybrid population, thus promoting a rapid accumulation of repetitive elements through subsequent segmental duplications or other recombination‐based mechanisms. Although speculative at present, centromere drive could commonly affect the evolution of admixed genomes, including the build‐up of reproductive isolation (e.g., Trier, Hermansen, Sætre, & Bailey, 2014) or the ancestry‐biased accumulation of repetitive elements (this study) during hybrid speciation.

Conclusions

Copy‐number variation in natural populations can contribute to the raw material for evolutionary change. The occurrence of high‐frequency gene duplicates with potentially adaptive gene dosage effects in invasive Cottus is in line with such a scenario but appears to be relatively rare in this study. We speculate that centromere drive followed by segmental duplications might contribute to the conspicuous accumulation of repetitive elements. While this appears to be attributable to hybridization and genomic admixture, it also suggests that it may be particularly pronounced in gene poor regions around the centromeres. On the one hand, this means that the massive copy number increase is unlikely to affect a large number of coding genes, while on the other hand, it makes it plausible that the large number of new repetitive element copies was not removed by selection. This warrants more studies that distinguish whether increased transposition rates or transposition‐independent mechanisms cause the accumulation of copy numbers in hybrid lineages and, more broadly, whether this may lead to rapid evolution of hybrid species.

Data Accessibility

Microarray raw data, including the aCGH calibration data, were archived at the Gene Expression Omnibus database under Accession no GSE93065 (http://www.ncbi.nlm.nih.gov/geo/). Raw sequence reads of Cottus genomes were archived at the European Nucleotide Archive (http://www.ebi.ac.uk/ena/) under Accession nos ERR1894831 to ERR1894860 (project Accession: PRJEB20064).

Author Contributions

S.D. and A.N.W. designed the experiments, analysed the data and wrote the manuscript. J.A. conducted the Illumina sequencing, and F.S. helped with all steps of the genome assembly and genome mapping. E.I. and X.L. contributed to the analysis of expression changes of genes and copy‐number changes of repetitive elements. All authors revised and approved the final manuscript.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.