Chromosomal Inversion Polymorphisms in Two Sympatric Ascidian Lineages.

Chromosomal rearrangements can reduce fitness of heterozygotes and can thereby prevent gene flow. Therefore, such rearrangements can play a role in local adaptation and speciation. In particular, inversions are considered to be a major potential cause for chromosomal speciation. There are two closely related, partially sympatric lineages of ascidians in the genus Ciona, which we call type-A and type-B animals in the present study. Although these invertebrate chordates are largely isolated reproductively, hybrids can be found in wild populations, suggesting incomplete prezygotic barriers. Although the genome of type-A animals has been decoded and widely used, the genome for type-B animals has not been decoded at the chromosomal level. In the present study, we sequenced the genomes of two type-B individuals from different sides of the English Channel (in the zone of sympatry with type-A individuals) and compared them at the chromosomal level with the type-A genome. Although the overall structures were well conserved between type A and type B, chromosomal alignments revealed many inversions differentiating these two types of Ciona; it is probable that the frequent inversions have contributed to separation between these two lineages. In addition, comparisons of the genomes between the two type-B individuals revealed that type B had high rates of inversion polymorphisms and nucleotide polymorphisms, and thus type B might be in the process of differentiation into multiple new types or species. Our results suggest an important role of inversions in chromosomal speciation of these broadcasting spawners.

Significance

Chromosomal rearrangements, especially inversions, are considered to play a major role in local adaptation and speciation in various taxa. In the present study, we sequenced the genomes of two individuals of a species of ascidian (an invertebrate chordate), and found high rates of inversion polymorphisms and nucleotide polymorphisms. In addition, alignments against chromosomes of a closely related ascidian revealed many inversions differentiating these two lineages of ascidians. Our study raises a possibility that inversions may have played an important role in chromosomal speciation in ascidians, and provides genomic resources for investigating chromosomal rearrangements in local adaptation and speciation.

Introduction

Chromosomal rearrangements, especially inversions, have an important role in the origin of new biological diversity (Sturtevant and Dobzhansky 1936; Dobzhansky and Pavlovsky 1958; Carson 1973; Rieseberg 2001; Wellenreuther and Bernatchez 2018). Such rearrangements locally suppress meiotic recombination in heterozygotes and thereby prevent gene-flow in the rearranged regions. Such gene-flow prevention may promote local adaptation, which may facilitate heterosis. In addition, chromosomal rearrangements may reduce the fitness of heterozygotes, and can thereby become a reproductive barrier. Inversions in particular have been studied in many organisms including land plants, insects, birds, and mammals (Carson 1973; Noor et al. 2001; Hoffmann and Rieseberg 2008; Lemaitre et al. 2009; Ayala et al. 2013; Twyford and Friedman 2015; Kupper et al. 2016; Tuttle et al. 2016).

Meanwhile, many aquatic animals obtain fertilization by passive dispersal of water-borne sperm. This reproductive mode differs from copulatory internal fertilization seen in insects, birds, and mammals in ways that might have important consequences for genetic diversity and speciation (Pogson 2016). Ascidians, which are marine invertebrate chordates, provide an ideal model system to study chromosomal evolution in broadcast-spawning aquatic animals, because genome sequences have been decoded in two species of the genus Ciona (C. intestinalis and C. savignyi) (Dehal et al. 2002; Vinson et al. 2005; Satou et al. 2008).

The species C. intestinalis contains at least two cryptic species, which have been identified by comparative genomics, mitogenomics and ecological studies (Suzuki et al. 2005; Kano et al. 2006; Caputi et al. 2007; Iannelli et al. 2007; Nydam and Harrison 2011). These two animal types, “type A” and “type B,” exhibit different susceptibility to temperature and morphological differences (Sato et al. 2012; Sato et al. 2015; Malfant et al. 2017). A recent study (Brunetti et al. 2015) has proposed to recognize them as distinct species and to call the former C. robusta and the latter C. intestinalis. However, to avoid possible confusion, since the first Ciona genotype to be published (type A) was referred to C. intestinalis (Dehal et al. 2002; Satou et al. 2008), we will use “type A” and “type B” in the following text.

These two types of Ciona occupy distinct geographical regions; type A lives mainly in Pacific and Mediterranean Sea, and type B lives only in the North Atlantic. A small sympatric area occurs in the western English Channel and southern Brittany, where natural hybrids of type-A and type-B animals are found at low frequency (Sato et al. 2012; Bouchemousse et al. 2016). We here provide chromosomal assemblies of two type-B individuals collected on the north and south sides of the English Channel, where type-A and type-B animals are in sympatry, to study chromosomal rearrangements in these animals.

Results and Discussion

Genome Size Estimation for Two Type-B Individuals

To obtain chromosomal-level assemblies for the two type-B individuals, sampled at Roscoff (France; specimen R) and Plymouth (England; specimen P), we first estimated their genome sizes to be 124 and 136 Mb, within which 95 and 91 Mb were respectively unique, using Illumina short-sequencing reads from sperm DNA of these two specimens by a k-mer profile method (supplementary table S1, Supplementary Material online). Thus, although these genome sizes are different, it is likely that this difference is mainly due to repeat sequences.

A previous report estimated that genomes of the type-A HT-line and two wild caught type-A individuals are 114–120 Mb in length (Satou et al. 2019). From the same data sets (DRA002218 and DRA008508) (Satou et al. 2015, 2019), we estimated their unique lengths to be 94–98 Mb. Therefore, it is likely that type-A and type-B individuals have nonrepeat regions of similar length.

Although the heterozygosity rate of type-A animals has been estimated to be 1.1–1.2% (Dehal et al. 2002; Satou et al. 2012), those of specimens R and P were 3.0% and 3.6%, respectively. Our data thus indicated that genomes of type-B animals were more heterozygous than those of type-A animals, consistent with transcriptome data indicating that exonic sequences of type-B animals are more polymorphic than those of type-A animals (Roux et al. 2013).

Constructing Assemblies for Two Type-B Individuals

To construct chromosomal assemblies for the two type-B individuals, we employed the strategy shown in supplementary figure 1A, Supplementary Material online. First, for each specimen, we built contigs from long Nanopore reads with the NECAT assembler (Chen et al. 2021) and polished obtained contigs with Illumina reads using the Nextpolish software (Hu et al. 2020). The total lengths of the contigs greatly exceeded the estimated genome size. Because the two specimens are diploid and alignments among contigs indeed showed many overlaps, we put such overlaps aside using the purge_dups software (Guan et al. 2020), and obtained nonoverlapping contigs. Both of these contig sets were 140 Mb in length (table 1). The contigs put aside were kept for studying haplotypes. We hereafter call the former contig sets the main contig sets, and the latter the haplotype contig sets. The two main contig sets were then aligned against the type-A chromosomes (Satou et al. 2019). Most contigs were aligned to unique positions of the type-A genome (supplementary fig. S2, Supplementary Material online).

Table 1

Basic Statistics of the Contigs and Scaffolds

		Specimen R		Specimen P
		Contigs	Scaffolds	Contigs	Scaffolds
Total nucleotide length (bp)		140,229,552	140,258,052	140,852,390	140,915,890
The number of “N”s in contigs/scaffolds (bp)		436	28,936	525	64,025
Number of contigs/scaffolds		92	35	234	107
Total nucleotide length in chromosomes		N.A.	134,327,883	N.A.	128,639,392
N50 (bp)		2,616,637	9,879,317	1,343,034	8,719,396
L50		19	6	34	7
N90 (bp)		921,931	6,498,296	342,918	5,508,368
L90		53	13	113	14
BUSCO score	Complete	92.8%	92.9%	93.6%	93.7%
	Fragmented	0.8%	0.9%	1.2%	1.1%
	Missing	6.4%	6.2%	5.2%	5.2%

The alignments of the longest contigs of specimens R and P against the type A genome are shown in figure 1. Although small inversions and one possible insertion were observed, the entire contigs were aligned with chromosomes 10 and 8, respectively. These observations indicate that genomic structures are globally conserved between type-A and type-B animals. Therefore, we constructed scaffolds from the alignments shown in supplementary figure S2, Supplementary Material online.

Comparisons among genomic sequences for two type-B animals and an inbred type-A animal. (A, B) Alignments of the longest contigs of specimens (A) R and (B) P against type-A chromosomes. An arrowhead indicates an insertion in specimen R, and arrows indicate inversions. Alignments of all contigs with type-A chromosomes are presented in supplementary figure S2, Supplementary Material online. (C–E) Exons, introns, 1-kb upstream regions, and intergenic regions including 1-kb upstream regions of the genome of the type-A HT line were mapped to the genomes of specimens R and P. Proportions of mapped nucleotides (C), nucleotide identities in mapped regions (D), and mean gap numbers per sequence (E) are shown.

N50, N90, and other related values were greatly improved in the resultant scaffold sets (table 1). The BUSCO program (Simao et al. 2015) indicated that over 93% of “metazoan genes” were included in the assemblies for specimens R and P. Notably, all these metazoan genes found by BUSCO were found in the chromosomes, and not in contigs that were not associated with the chromosomes.

We also confirmed that the animals we used were type B by molecular phylogenetic analyses using five loci that have been reported to be divergent between two types (Nydam and Harrison 2011) (alignments are shown in supplementary figs. S3 and S4, Supplementary Material online).

We also mapped 318 genes for transcription factors and signaling molecules of the type-A genome to the genomes of specimens R and P, because these genes were well annotated (supplementary fig. S5, Supplementary Material online). As was expected from the genomic alignments, chromosomal positions of regulatory genes were largely conserved, although two possible inversions and a possible interchromosomal translocation were found (see below for analyses of such structural variations). Among these genes, 317 and 311 genes were successfully mapped to chromosomes of specimens R and P, respectively. Thus, >98% of genes are present in the assembled chromosomes of these specimens.

Comparisons among the Genomes of the Type-A and Type-B Animals

To compare the genomes of the type-B animals and the type-A HT-line at the nucleotide level, we split the entire genome sequences into 500-bp-long fragments, and aligned them to the genome sequences of different animals by the BLAT program (Kent 2002) with default parameters (supplementary fig. S6, Supplementary Material online). These fragments were more alignable between the type-B animals than between the type-A and type-B animals, and nucleotide identities within aligned fragments were higher between the type-B animals than between the type-A and type-B animals. These observations indicate that the type-B genomes are more similar to each other than to the type-A genome.

We next used a gene model set that we made previously for the type-A genome (Satou et al. 2019), and mapped exons, introns, 1-kb upstream regions, and intergenic regions (including 1-kb upstream regions) to the type-B genomes by BLAT (Kent 2002). Proportions of nucleotides successfully aligned indicate that exons are more conserved between the type-A and type-B genomes than introns and intergenic regions (fig. 1). Although identities of the aligned nucleotides were almost the same among the exons, introns, upstream regions, and intergenic regions (fig. 1), numbers of gaps found in the alignments showed large differences (fig. 1). Thus, exons are highly conserved between type-A and type-B animals, whereas intergenic regions are highly varied between the two types.

Occurrence of Chromosomal Inversions in the Type B Genomes

To examine structural variations, we mapped the KY gene model set of the type-A genome to the genomes of specimens R and P. Among 14,072 gene models, 8,479 and 8,524 models were mapped with high scores onto the type-B genomes of specimens R and P, respectively. The order of genes on these three genomes was highly conserved whereas structural variations were also detected (supplementary fig. S7, Supplementary Material online). This contrasts with extensive intrachromosomal rearrangements between the type-A animals and C. savignyi (Hill et al. 2008; Satou et al. 2019).

To identify translocations, we searched the genomic sequences for blocks containing three or more genes that were not found in the positions expected from the type-A genome. We found two small interchromosomal translocations between the genomes of specimen P and type A, which were also confirmed with the genomic alignment (supplementary figs. S7–S9 and table S2, Supplementary Material online). On the other hand, we did not find such translocations between the genomes of specimen R and type A.

Similarly, to identify inversions, we searched the genomic sequences for blocks containing three or more genes that were mapped in the reverse direction to the order in the type-A genome. We found 21 and 20 inversions in the genomes of specimens R and P, respectively (supplementary fig. S7, Supplementary Material online). Among them, 15 sites were common, and the remaining 11 sites were specific to specimen R or P (supplementary fig. S10 and tables S3 and S4, Supplementary Material online). Thus, our data indicate that there are structural variations not only between type-A and type-B animals but also between specimens R and P.

Structural Variations among the Genomes of Type B Animals

To experimentally verify structural variations between specimens R and P, we focused on an inversion in chromosome 7 because it contained the largest number of genes over an ∼900-kb region (fig. 2). We found R-derived haplotype contigs that were structurally different from the chromosome of specimen R, but the same as both the specimen-P and the type-A chromosomes. This observation indicated a possibility that the genome of specimen R is heterozygous in this region. Genomic PCR using four primers flanking the junctional sites (fig. 2) demonstrated that specimen R was indeed heterozygous in this region and that specimen P was homozygous (fig. 2).

Experimental verification of three inversions. Genomic alignments of regions of chromosomes (A) 7, (C) 3, and (E) 1 of specimen R against corresponding chromosomal regions of specimen P, corresponding haplotype contigs of specimens R and P, and corresponding type-A chromosomes. Forward-to-forward alignments are shown by magenta, and forward-to-reverse alignments are shown by cyan. Illustrations below the alignments indicate PCR primer positions for experimental verification of the genomic arrangement of these regions. These primers were used to examine which set gave specific amplification. PCR products for chromosomes (B) 7, (D) 3, and (F) 1 were analyzed by agarose gel electrophoresis.

We also confirmed experimentally an inversion on chromosome 3, to which seven genes were mapped. The genomic alignments indicated that specimen P was heterozygous whereas both haplotypes of specimen R shared the same structure with type A (fig. 2). A PCR analysis demonstrated that specimen P was indeed heterozygous in this region (fig. 2).

Because small variations or variations in genomic regions where genes are rarely encoded or less-conserved genes are encoded may not be detected by the above method, we manually inspected genomic alignments of the two type-B specimens against the type-A genome (supplementary fig. S2, Supplementary Material online), and found nine additional inversions (supplementary table S5, Supplementary Material online). Four and two inversions were specific to specimen R and specimen P, respectively. We experimentally verified one of the newly identified inversions on chromosome 1 (∼250 kb) (fig. 2). That is, specimen R was heterozygous, and specimen P was homozygous for this region (inverted against the type-A genome).

To understand how prevalent the inversions we found were in type-B populations, we collected eight wild specimens at Roscoff and Plymouth, RO1–RO4 and PL1–PL4, respectively, and performed genomic PCR for the same three inversions (fig. 3). For the inversion sites in chromosomes 3 and 7, all samples were homozygous and all haplotypes were in the same direction (fig. 3). The directions of these two sites of all eight samples were the same as those of type A. For the inversion site in chromosome 1, we found one heterozygous individual, which was collected at Plymouth (fig. 3). The direction of this site in six of the other seven homozygous specimens was opposite to that of type A. Thus, type-B animals indeed have inversion polymorphisms. It should be highlighted that the number of inversions in type-B populations is likely to be underestimated in the current analyses, because our analyses only used gene orders and genomic alignments 3 kb or more in length and because our analyses were limited to two type-B specimens.

Inversions in wild animals. (A–F) Genomic PCRs for three genomic positions on chromosomes (A, D) 7, (B, E) 3, and (C, F) 1 using DNA from animals collected at (A–C) Roscoff and (D–F) Plymouth. Primers were the same as the ones used in figure 4, and primer pairs are shown below the photographs. As controls, genomic PCRs using DNA from specimens R and P are also included. (G–I) Genomic PCR results for the above three positions are respectively illustrated with the arrangement in the type-A genome.

Note: Figure Replacement Requested.

On the other hand, it is not likely that type-A animal populations extensively retain inversion polymorphisms, for the following two reasons. First, only one inversion (in chromosome 4) of the HT-line was identified by a comparison between the HT-line genome and the KH version of the genome that is derived from a different type-A individual (Satou et al. 2019). Second, mapping of Hi-C assay data, which are derived from wild animals, to the HT-line genome identified only one inversion at the same location in chromosome 4 (Satou et al. 2019). Therefore, it is likely that type-B animals have more inversion polymorphisms than type-A animals.

Inversions may reduce fitness of hybrids postzygotically and reduce efficiency of meiotic recombination in heterozygotes (Zetka and Rose 1992; Dresser et al. 1994; Lamb et al. 2007; Massip et al. 2010; Kirkpatrick et al. 2012; Ederveen et al. 2015). If reduced fitness of hybrids indeed contributes to speciation between sympatric type-A and type-B animals, the high inversion polymorphism rate in type-B populations may indicate that this group may be separating into further subtypes. Indeed, the high heterozygosity rates documented by the present study and previous studies (Tsagkogeorga et al. 2012; Roux et al. 2013) indicate high nucleotide diversity in type-B animals, which supports this idea.

Inversions locally suppress meiotic recombination in heterozygotes. The inversions we found in the present study are relatively small (3 kb ∼ 873 k; supplementary tables S3 and S5, Supplementary Material online). Because the recombination rate in Ciona is estimated to be 25–49 kb/cM (Kano et al. 2006), these inversions correspond to 0.06–35 cM. Although double crossover may yield high levels of gene flux (allele exchange during meiosis in heterokaryotypic oogenesis) in central regions of long inversions (Navarro et al. 1997), the inversions we found are relatively small and therefore high levels of gene flux are not expected. As a result, inversions may play a role in preventing gene flow in inverted regions between type-A and type-B animals, and this prevention may promote local adaptation and speciation (Noor et al. 2001; Rieseberg 2001; Wellenreuther and Bernatchez 2018). Namely, it is possible that an inverted region contains genes involved in local adaptation or speciation and that these genes are protected from gene flows between type-A and type-B animals. Candidate genes, which are contained in the inverted regions, are listed in supplementary tables S4 and S5, Supplementary Material online. Because intergenic regions are more highly variable than coding regions (fig. 1), it is possible that changes not only in coding regions but also in regulatory elements may also be involved.

Although the significance of inversions has been documented in insects (Sturtevant and Dobzhansky 1936; Dobzhansky and Pavlovsky 1958; Carson et al. 1967; Carson 1973; Noor et al. 2001; Ayala et al. 2013), birds (Kupper et al. 2016; Tuttle et al. 2016), and mammals (Lemaitre et al. 2009), our data highlight the likely importance of inversions in speciation of Ciona, an invertebrate with mating governed by interactions of aquatic gametes. In addition, our genomic data provide a platform to dissect molecular bases for speciation at the nucleotide sequence level. Indeed, structural variations we found includes small variations that are probably impossible to find by cytological analyses.

In the apparent absence of extrinsic prezygotic and postzygotic barriers to hybridization, there is a very low rate of contemporary hybridization between types A and B in sympatry (Sato et al. 2012; Malfant et al. 2017). The explanation for this has been sought in various suggested intrinsic postzygotic barriers, including genomic incompatibilities in second-generation hybrids (Dobzhansky–Muller incompatibilities) (Roux et al. 2013; Bouchemousse et al. 2016; Malfant et al. 2018) and mitochondrial-nuclear incompatibility (Ohta et al. 2020). We suggest that the influence of linkage groups associated with frequent chromosomal inversions is an additional contributing factor, as was argued for chromosomal evolution of Ciona species based on comparison of the genomes of C. intestinalis type-A (C. robusta) and C. savignyi (Satou et al. 2019).

Materials and Methods

Biological Materials

Two C. intestinalis type-B individuals were obtained from wild populations near Roscoff and Plymouth marine stations, respectively, for genome sequencing. Additional wild animals were also collected at these sites.

To confirm that the animals we used for genomic sequencing were type B, we used genomic sequences of five loci, Fgf4/5/6 (this gene annotation was likely incorrect, because the sequences found in the public database were all mapped to a region within an intron of the Fgf receptor gene), Foxa.a (fkh), Jade, Patched, and Vesicular acetylcholine transporter (vAChTP); these loci have been reported to be divergent between the two types (Nydam and Harrison 2011). Sequences retrieved from the public nucleotide database were aligned using the Clustal Omega program (Sievers et al. 2011), and then alignments were manually adjusted. After removing gaps, alignments were used to construct molecular phylogenetic trees by the maximum likelihood method with the PhyML program (Guindon and Gascuel 2003). Trees were tested with 100 bootstrap pseudoreplicates.

Genome Sequencing

In Nanopore sequencing for specimen R, a 1D Ligation Sequencing Kit, SQK-LSK109 (Oxford Nanopore Technologies) was used to prepare a library. The genome DNA was sheared by Megarupter (Diagenode) and sequenced on a PromethION with a R9.4.1 flowcell. In Nanopore sequencing for specimen P, a Rapid Sequencing Kit (SQK-RAD004, Oxford Nanopore Technologies) was used to prepare a library, and the library was sequenced on a MinION with a R9.4.1 flowcell. For Illumina sequencing, genome DNA was extracted from specimens R and P using phenol chloroform and InnuPure C16 touch with smart DNA prep Kit (Analytik Jena). We used NEBNext Ultra II FS DNA Library Prep Kit for Illumina (New ENGLAND BioLabs) and KAPA HyperPlus Library Preparation Kit (Roche). The Illumina platform (Hiseq2500 and NovaSeq6000) was used for sequencing. PacBio RSII sequencing involved an SMRTbell Template Prep Kit 1.0 (Pacific Biosciences) to prepare a library, which was sequenced using a PacBio RSII sequencer employing P6-C4 chemistry with 360-min movie lengths.

Contig assembly was performed using the NECAT program with default parameters (Chen et al. 2021). Then contigs were polished using Nextpolish (Hu et al. 2020) with Illumina reads obtained from the same specimens. Possible haplotypic contigs were identified using the purge_dups program (Guan et al. 2020); for specimen R, PacBio reads were used, and for specimen P, the aforementioned Illumina reads were used. Contigs were aligned with the genome of an inbred type-A animal (HT-line) (Satou et al. 2019) using Nucmer (Kurtz et al. 2004). After manual inspection, contigs were joined. We inserted 500 “N”s between contigs.

Genome Size Estimation

Genome sizes of specimens R and P were estimated by counting all possible 21-mers with Jellyfish (Marcais and Kingsford 2011) and analyzing them with Genomescope (Vurture et al. 2017).

Genomic Alignments

Pairwise alignments between the type-A genome and the genome of specimen R or P were performed with Nucmer (Kurtz et al. 2004). Short alignments of <1,000 bases were removed. The entire coding regions, individual exons, introns, upstream regions, and intergenic regions were deduced from the KY gene model set, which were made for the type-A genome (Satou et al. 2019), and these nucleotide sequences were mapped with BLAT (Kent 2002). The order of genes in the type-A genome was compared with the order of genes mapped to the type-B genomes.

Experimental Validation of Inversions

To confirm inversions found in the genomic assemblies, we performed PCR. Primer sequences are shown in supplementary table S7, Supplementary Material online. PCR was performed with Primestar Gxl (Takara Bio).

Supplementary Material

Supplementary data are available at Genome Biology and Evolution online.

Click here for additional data file.