Phylogeographic study using autonomous reef monitoring structures indicates fast range expansion of the invasive bryozoan Juxtacribrilina mutabilis

This is a phylogeographical study of Juxtacribrilina mutabilis, a recently described bryozoan from Japan with sightings in Sweden, Norway, and Maine (US), to test how fast and far it has expanded across Europe in recent years. J. mutabilis settles easily on ship hulls, making it a useful model organism for studying long-distance invasion pathways. The study was conducted using Autonomous Reef Monitoring Structures (ARMS) to monitor the entire European coastline and Svalbard for J. mutabilis using DNA metabarcoding. During the time between its first sighting in Europe 2008 and the last ARMS retrieval 2020, the data shows how J. mutabilis has gained a pan-European distribution, being genetically identified in 14 new distinct locations. Presence/absence of barcodes were confirmed by image-based identification in 74% of the samples. Fourteen haplotypes never reported before were discovered in a 264 bp region of the cytochrome oxidase I gene. Two haplotypes (HP1 and HP3) occurred frequently and geographically widely dispersed, indicating intercontinental connectivity. Two locations, Koster and Getxo showed particularly high genetic diversity with similar haplotype networks suggesting continuous gene flow across oceanographically unlinked regions. Given the recent description of J. mutabilis and the relatively few historical encounters, the genetic diversity described here suggests an unusually fast range expansion within the last two decades. Such global spreading events of fouling organisms may become more common in the future as a result of more frequent and interconnected ship traffic.


Introduction
Juxtacribrilina mutabilis (Ito et al., 2015), previously Cribrilina mutabilis, is an invasive cheilostome bryozoan.It was first described at lake Akkeshi (Japan) and has also been found in Kristineberg (Sweden), Bergen (Norway), and in Cosco Bay, Maine (US) (Ito et al., 2015).Before that it is thought to have been observed under the name Membraniporella aragoi (Audouin, 1826) in the Russian Northwest Pacific by Audouin (1826) and Kubanin (1975Kubanin ( , 1977Kubanin ( , 1997) ) (Dick et al., 2020).The north-eastern passage is a possible vector of invasion from Asia to Scandinavia, from where it most Vol:.( 1234567890) likely has spread to Maine.Its main reported impact is the fouling of eelgrass, which is its most common settling substrate, but it has also been found on other Zostera species, several algae, plastic panels and strips and ship hulls (Dick et al., 2020).
The ability to settle on ship hulls in conjunction with previous demonstrated spreading over large distances makes J. mutabilis a good model for studying invasion pathways and expansion rates associated with new shipping routes (Gunnarsson, 2021).In addition to new shipping routes, also the intensity of ship traffic has multiplied over the past decades (UNCTAD, 2018).Hence, we can expect to see new patters of spreading, not only with regard to the number of species but also in the rate of expansion.Estimates of expansion rate can have important applications in conservation management.Quantitative measurements of invasion speed are difficult to produce as they are likely dependent on underlying evolutionary processes (Philips, 2015).Yet, understanding the scale at which species can invade regions will be important for defining windows of opportunity for counteraction.With climate change shown to cause profound marine species distribution shifts (Perry et al., 2005) we can expect even more unwarranted invasions in the near future, especially alongside ship traffic that connects previously incompatible ecosystems.
In this study we investigate the phylogeographic structure of J. mutabilis to test how far and how fast the species has expanded in European waters in the past decade.More specifically, we describe the species current distribution and genetic diversity in Europe, try to identify signatures of recent translocation of the species across the Pacific and Atlantic Ocean, and investigate its true place of origin.To this end we utilized the genetic monitoring data collected the European ARMS MBON program (www.armsmbon.eu/), which deploys settling plates, so-called Autonomous Reef Monitoring Structures (ARMS) in ports, marinas, and nature reserves along the European coastline and in the polar regions.These plates are used to collect eDNA metabarcoding data (Ficetola et al., 2008) as well as photographic images of the settlement plates and allow for analysis of genetic diversity across benthic habitats in European waters (Obst et al., 2020).

Analysis of ARMS samples
Sixty-six ARMS samples for which all data are openly accessible were analyzed across 18 locations in Europe (Table 1).Samples were deployed, retrieved, photographed, and processed by individual partners in the ARMS-MBON network (www.armsmbon.eu), while DNA extraction, PCR amplification and sequencing was done centrally following the Molecular Standard Operating Procedures (MSOP) published by the ARMS-MBON program (Obst et al., 2020), using the mlCOIintF and jgHCO2198 primers.To confirm any genetic match, the plate photographs of all 66 ARMS were scanned for presence/absence of J. mutabilis.In some cases, J. mutabilis was not identified on the images but on closer inspection there were instances where ancestrulas (initial stage of bryozoan colony development) were detected on the images.Ancestrulas could not be determined to species level, yet could potentially belong to J. mutabilis.The number of colonies per sample was not taken into account, as deployment time varied between the samples (Fig. 1).
Bioinformatic analysis was performed with a custom pipeline written in the R programming language (R-Core-Team, 2019) available in the Supplementary material.Raw sequencing files in fastq format (ENA accession numbers are given in Table 1) were processed with Dada2 (Callahan et al., 2016).Primer location and orientation was assessed, while primer removal was performed with Cutadapt (Martin, 2011).Quality profiles for all samples and markers were assessed with FastQC.Low-quality read ends (quality score < 2) were trimmed off.Furthermore, the error model was trained and evaluated, resulting in a good match between observed and expected errors.Thereafter de-replication and inference of Amplicon Sequence Variants (ASV) were performed, merging forward and reverse reads.Finally, chimeras and singletons were removed using "pseudo-pooling."Remaining sequences were trimmed to a minimum length to exclude short reads (< 300 bp).After filtering, the number of output reads per sample varied widely from none to tens of thousands.Samples with less than 5,000 reads were considered to be of "poor quality," but remained in the analysis.Taxonomic assignment of COI ASVs was done against the BOLD reference library (Ratnasingham & Herbert, 2007) using BOLDigger (Buchner & Leese, 2020).Matches of J. mutabilis of 97% or higher were singled out from the data and tracked throughout the various locations.The 97% threshold was chosen empirically, as most matches were > 99% or < 90%.Each ASV was checked for alternative matches manually in the BOLD database to rule out any incorrect taxonomic assignments.ASVs were screened for possible pseudogene copies using MACSE v2.05 (Ranwez et al., 2018) and the 'enrichAlignment' function with genetic code 5, allowing no in-frame stop codons, frameshifts, and a maximum of three deletions on the amino acid level.The reference alignment consists of curated COI reference sequences mined from BOLD, representing the sequence diversity of deposited bryozoan COI sequences-and is provided by the developers of MACSE (Delsuc & Ranwez, 2020).All ASVs could be aligned to the reference alignment, and hence no ASVs were considered pseudogenes.All ASVs matching J. mutabilis were submitted to the ASV portal of the Swedish Biodiversity Data Infrastructure (asv-portal.biodiversitydata.se/) and are also made available to GBIF (Martaeng et al., 2022).

Haplotype network analysis and haplotype map
In order to analyze frequency and relationships between the discovered genetic sequences, a haplotype network was calculated combining the public J. mutabilis sequences (downloaded from BOLD) with all unique COI sequences found in the ARMS data set.The method was the same as in Ito et al. (2015), a Tight Span Walker network created using PopArt (Leigh & Bryant, 2015), downloaded at http:// popart.otago.ac.nz/ downl oads.shtml (last accessed December 2021).The original study used a region of 555 bp length (some 631 bp long), while the ARMS data only covered a 313 bp region.The overlap between the public sequences and the ARMS sequences was 264 bp.Hence, all sequences were aligned and trimmed in R and a haplotype network was calculated for the overlapping 264 bp region.In order to compare the haplotype composition between the two locations with highest diversity, two additional networks were calculated for the two locations with the highest genetic diversity, Getxo and Koster using the 264 bp alignment.Finally, haplotype diversity was plotted on a map using QGIS (QGIS Development Team, 2022).

Statistical analysis
Statistical analysis was made for the whole dataset as well as for Getxo (ES) and Koster (SE) individually due to their notable haplotype diversities.First, measurements of DNA polymorphism; haplotype diversity Hd and nucleotide diversity Π (Nei, 1987), were calculated.Next, the following statistical tests were calculated: Tajima's D-test (Tajima, 1989), Fu's Fs test (Fu & Li, 1993), Ramos-Onsins, and Rozas's R2 test and raggedness index (rg) (Ramos-Onsins & Rozas, 2002).All statistical calculations were done in DnaSP v6.12.03 (Librado & Rozas, 2009), utilizing coalescent simulations of 1000 re-sampling replicates.The null hypothesis for the tests was that the population is stable, while the alternative hypothesis is that the population is (rapidly) expanding.

Results
Of the 66 ARMS samples investigated, 9 were of poor genetic quality (delivering < 5,000 sequence reads) but could still be positive for J. mutabilis (Table 1).The species was identified genetically and photographically in 23 samples, while we found only genetic signals in 13 samples, and only photographic evidence in 4 samples, resulting in overall 40 positive samples.The remaining 26 samples had no matches at all.Photographic presence/absence and genetic presence/absence was in concurrence in 49 of the 66 samples (74%).Overall, J. mutabilis was identified genetically in 14 new locations showing a pan-European distribution from Svalbard in the North to the Eastern Mediterranean in the South, and from the Atlantic coast in the West to the Gulf of Finland in the East (Fig. 2).
Overall, we retrieved 21 unique ASVs with a length of 313 bp from all ARMS samples.When aligned with the public sequences, the resulting 264 bp fragment returned 16 haplotypes (Fig. 2).In this data set, three haplotypes were previously known (HP1, HP2, and HP3), but only two of the sequences were recovered in the ARMS data (HP1, HP3).One (HP2) was not found in our data set.Hence, we included 17 globally known haplotypes in our analysis of which 14 were newly discovered by this study (HP4-HP17).The complete haplotype network (Fig. 2) further expands on the network produced by Dick et al. (2020) with emphasis on European coastal water.The genetic centroids of the network (HP1 and HP3) are also the most frequent and geographically widespread haplotypes (Fig. 2; Table 2).HP1 is present at six locations while HP3 is present at every location except Akkeshi (JP).
Two locations, Getxo (ES) and Koster (SE) have notably high haplotype diversities: 7 and 8 haplotypes respectively (Table 2).The individual haplotype networks calculated for these locations (Fig. 3) are very similar, despite their geographical distance.Both networks share the same core structure with a few variant haplotypes in the periphery.The statistical analysis of the DNA polymorphism (Table 3) showed high haplotype diversity but low nucleotide diversity for all three sets (entire geographic range, Getxo, and Koster).All tests were, however, statistically insignificant in relation to various thresholds commonly associated with these tests and hence failed to reject the null hypotheses that the population is stable.

Biogeography and range expansion
In this study we mapped distribution and genetic diversity of the recently described invasive species from Japan, Juxtacribrilina mutabilis, across Europe using standardized genetic sampling methods.Our results show considerable genetic diversity of J. mutabilis across its entire European distribution range.The high haplotype diversity in conjunction with low nucleotide diversity observed at Koster, Getxo and in all samples combined (Table 3) suggest a rapid recent expansion, i.e., a "population bottleneck followed by rapid population growth and accumulation of mutations," according to Grant and Bowen (1998).But the statistical tests (e.g., Tajima's D) failed to reject the null hypothesis that the population(s) are stable.A reason of the weak performance of the statistical tests may be the relatively short fragment size and the small sample sizes, especially for the locations Koster and Getxo.
The similar genetic signatures across distant locations described in this study (Figs. 2, 3; Tables 2, 3) Table 2 Locations for each COI haplotype found in the dataset and in the original tissue samples Only locations with at least one haplotype are shown.Names are abbreviated: Akkeshi (Akk), Angbat (Ang), Casco Bay (Cas), Crete (Cre), Getxo (Get), Gothenburg (Got), Helsingborg (Hel), Knokke-Heist (Kno), Koster (Kos), Kristineberg (Kri), Kummelkobben (Ku), Laesoe (Lae), Lysekil (Lys), Marstrand (Ma), Plymouth (Ply), Roscoff (Ros), Svalbard (Sva) and Toralla (Tor).Also visualized in a haplotype map, Fig. 2 HP Akk Ang Cas Cre Get Got Hel Kno Kos Kri Ku Lae Lys Ma Ply Ros Sva Tor are indicative of recent and intense translocations of the species.HP3 is present at every location except Akkeshi, while the previous genetic distance of 6 bp per 555 bp between Japanese and European haplotypes shown by Dick et al. (2020) has been reduced to no distance at all.HP1 has been found in Akkeshi, Getxo, and at several Swedish locations.This intercontinental distribution of haplotypes HP1 and HP3 suggests swift migration between Japan, Europe, and North America.Furthermore, the similar haplotype networks in Getxo and Koster indicate continuous gene flow across oceanographically unlinked regions in Europe, such as the Bay of Biscay and the Skagerrak.In both cases intense and regular ship traffic may provide a plausible explanation for the observed gene flow.
The species was first found in Europe in Bergen in 2008, but was not reported from Sweden until 2011 (Dick et al., 2020).Reasons for such non-observance could be taxonomic confusion with similar species such as Callopora rylandi or Cribrilina annulata caused by the polymorphism exhibited by J. mutabilis (Dick et al., 2021).Alternatively, J. mutabilis could have also been rare in Europe for a long time  and only recently became a native invasive.However, given the strong taxonomic expertise in the region we consider these scenarios unlikely and we see strong evidence that the species has invaded all regional seas of Europe (except the Black sea) within the 12 years between its first European sighting 2008 and the last retrieved ARMS data 2020.This rate of continental expansion outperforms several other marine invertebrate groups.Carcinus maenas for example showed its fastest regional expansion when it emigrated from San Francisco (California) to British Columbia over one decade (Carlton & Cohen, 2003).Neogobius melanostomus, another renowned invasive species, needed one decade to spread through the Great Lakes of North America, and approximately two decades to spread from the Baltic sea to the English Channel (Kornis et al., 2012).A third example is the invasive polychaete Marenzelleria viridis which was first sighted in Scotland in the late 1970s and which reached the Gulf of Bothnia and Gulf of Finland roughly two decades later (Maximov, 2011).All three species have a much shorter expansion range per decade than J. mutabilis.
In contrast to the recent expansion of J. mutabilis in European waters, the species may have been present in the Northwest Pacific for much longer, if the discoveries of Audouin (1826) and Kubanin (1975Kubanin ( , 1977Kubanin ( , 1997) ) are to be interpreted as J. mutabilis.Together, these historical records and the type locality in northern Japan suggest that the Northwest Pacific remains the most likely place of origin.However, in order be certain about the origin of the species, genetic diversity of J. mutabilis should also be assessed in the Northwest Pacific.
In conclusion, is appears that J. mutabilis has in less than two decades expanded from a presumed native range around the coast of Japan to a globally distributed species showing the capacity to tolerate a great range of water temperatures and salinities and traveling across large distances with the help of ship traffic.The presented rate of expansion exceeds that of most other benthic invertebrates and may be indicative for the increased global connectivity due to intensified ship traffic and opening of new trading routes in the arctic.We therefore argue that biological monitoring programs as well as mitigation policies need to prepare for faster biological invasions.
Vol.: (0123456789) metadata is likely part of several studies currently produced or reviewed.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made.The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material.If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.To view a copy of this licence, visit http:// creat iveco mmons.org/ licen ses/ by/4.0/.

Fig. 2 (
Fig. 2 (Top) Tight Span Walker haplotype network for Juxtacribrilina mutabilis using all 264 bp haplotypes found in our European eDNA metabarcoding dataset and GenBank.The sizes of circles are scaled to the number of samples.(Bottom) Haplotype map created in QGIS visualizing the locations of all known 264 bp haplotypes of J. mutabilis found in European

Fig. 3 a
Fig. 3 a Tight Span Walker haplotype network forJuxtacribrilina mutabilis using only 264 bp sequences found in ARMS from Koster (Sweden).b The respective haplotype network for ARMS from Getxo (Spain).In both networks, the sizes of circles are scaled to the number of samples, and haplotypes have been color coded to match Fig.2

Table 1
Genetic and photographic matches of Juxtacribrilina mutabilis Sample period is the deployment and retrieval of ARMS.ENA accession numbers are used by ARMS-MBON for recording sequencing runs, one number per sample and fraction size.Positive quality means that the sample generate at least 5,000 genetic reads after processing.Positive genetic match means that at least one of the reads are interpreted as J. mutabilis.Positive photo match means that J. mutabilis was observed at least once on photos of that ARMS Vol:.(1234567890)