Monthly eDNA Monitoring of an Invasive Bryozoan, Bugulina californica, in Seawater Using Species-Specific Markers

Simple Summary Bugulina californica, marine invasive bryozoan, is hard to monitor the biomass and presence because of their habitat in underwater. Additionally, they have life stage difficult to find such as larva, and we need an effective survey method to detect whole life stages for monitoring early invasion stage. Therefore, we tried to applied environmental DNA to monitor the monthly changes of B. californica in harbors of Korea. We collect seawater environmental samples and developed a molecular target species detection method to detect B. californica DNA of monthly changes. We analyzed the environmental samples using our molecular markers and calculated the DNA copies. We determined method of environmental DNA assay as effectiveness survey technique for marine invasive species which has a non-visual life stage and spatial changes of whole biomass. Abstract Environmental DNA (eDNA) method used by many ecologists as effective investigation tool can detect endangered species, rare species, and invasive species. In case of invasive species, eDNA method help to monitor the target species when the species was hard to detect through the traditional survey such as the early stage of invasion, low abundance, and larva or juvenile stage. The bryozoan, Bugulina californica, was known as a marine fouling invasive species in Korea since its first reported in 1978. This species expanded nationwide, and damages to ascidian aquaculture through attached on the ship hulls and artificial facilities. To monitor the distribution and biomass of invasive bryozoan, B. californica, the qPCR analysis of environmental DNA was performed on seawater samples from 12 harbors. In this study, we designed species-specific markers which can calculate the detected DNA copies of B. californica, and the presence and monitoring of this species can be more accurately estimated by environmental DNA analysis than by traditional survey, in which it is difficult to identify the species. Real-time PCR analysis using environmental DNA is an effective monitoring method that can determine both the distribution and the monthly change in biomass of B. californica in Korea.


Introduction
Ecologists require effective investigative tools and methods specifically designed to detect endangered species, rare species, and invasive species with low cost, minimum stress for the environment, and high potential for detection [1][2][3]. Traditional methodssuch as trapping, netting, and electrofishing in field surveys-have been used routinely for the detection of various species [4]. However, these methods are unsuitable for the detection of many rare species, especially in the aquatic environment, where organisms are hard to detected visually [5][6][7]. Therefore, molecular DNA-based tools applicable for

Field Survey and Specimen Morphological Examination
The field geological coordinates for this study were provided in Table 1. To improve the detection rates of B. californica, we conducted the field survey in summer (August, 2017) when their abundance observed in maximal condition. A total of 12 survey sites were selected based on the enough artificial substrates that various fouling organisms can attached, regular intervals, and covered three coastal.
To obtain the target species sample for molecular analysis and identification, B. californica was sampled by hand in the investigation sites (harbors) this species shown most abundant. B. californica samples attached on the artificial substrates, ropes, buoy, and concrete structures, were collected in February 2017 at depths ranging 0.5-5.0 m from the Gonghyeonjin (38 • 45.04 E) harbors ( Figure 1). The collected samples were immediately fixed with 95% ethanol, and stored at room temperature (25 • C). All of the colonies were examined and identified using scanning electron microscopy (SEM), and our sample images were analyzed [35][36][37]. Remaining samples were used in DNA analysis for development and validation of species-specific markers.

Molecular and Morphological Identification of B. californica 2.2.1. DNA Extraction, Amplification, and Sequencing
Total genomic DNA was isolated from each sample of three B. californica sample using a DNeasy blood and tissue DNA isolation kit (QIAGEN, Hilden, Germany), in accordance with the manufacturer's instructions. The quality and concentration of the extracted genomic DNA were determined using a Nanodrop ND-1000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA). All genomic DNA samples were stored at −20 • C until use.
The barcode region of the mitochondrial cytochrome c oxidase I (COI, 658 bp) gene was amplified and sequenced. The primers used for COI amplification were LCO1490 (forward) and HCO2198 (reverse) ( Table 2) [38]. PCR was conducted using a 25 µL reaction mixture, containing 2.5 µL of 10× Ex Taq Buffer containing 20 mM MgCl 2 (Clontech, CA, USA), 1 µL of 2.5 mM dNTPs (Clontech, CA, USA), 1 µL of each primer at 10 pmol, 1.5 µL of 150-250 ng/µL template DNA, 0.3 µL of 5 U/µL Taq polymerase (Clontech, CA, USA), and 17.7 µL of DEPC-treated water. The following PCR conditions were used: initial denaturation at 94 • C for 5 min; 35 cycles of denaturation at 94 • C for 30 s, annealing at 50 • C for 90 s, extension at 72 • C for 90 s; and final extension at 72 • C for 7 min. The PCR products were separated by electrophoresis on 1.5% agarose gel stained with ethidium bromide (Bioneer, Daejeon, Korea), and 2 µL of 100 bp DNA ladder (Elpis Biotech, Daejeon, Korea) was loaded to show the size of the proteins. The PCR products were directly sequenced in both directions with the primers used for amplification (Cosmogenetech, Seoul, Korea). The results were analyzed by BioEdit [39], to assess the quality of sequencing. Three COI sequences of B. californica from Gonghyeonjin, Ulsan and Dadaepo were registered in GenBank of NCBI (http://www.ncbi.nlm.nih.gov accessed on 22 May 2021), and the accession numbers were MZ209214, MZ209213, and MZ209215, in respectively.

Molecular Identification
To verify the selection of mitochondrial COI as appropriate gene region for designing the species-specific markers, we conducted the phylogenetic analysis using selected region. The sequences obtained from three B. californica samples were manually aligned using Clustal X [40]. Furthermore, we investigated the genetic distances and phylogenetic relationships of three Korean B. californica with voucher specimens of B. californica collected in Tongyeong (34 • 50 23.05 N, 128 • 25 12.58 E) registered in the Marine Bryozoans Resources Bank of Korea (resources number: MBRBK-177; GenBank accession number: MZ217204) and four other species accepted as the Bugulina genus, registered in the NCBI: B. fulva (KC129719 in UK), Bugula flabellata (AY633484; AY061749 in Nambia), B. turbinata (AY633481; AY633480; AY633482 in Canada), and B. simplex (AY633478 in USA). After all sequences were aligned, we identified a 507 bp region of COI in which all nucleotides overlapped. Genetic distances were calculated according to the Kimura 2-parameter (K2P) model [41], and bootstrap analysis was conducted with 1000 replicates. Phylogenetic analysis was performed using the neighbor-joining method with MEGA7.0 [42]. Terebratalia transversa (FJ196085 in USA), of the family Terebrataliidae, of order the Terebratulida, class Rhynchonellata, and the phylum Brachiopoda, served as outgroup for phylogenetic analyses. We designed species-specific primers amplifying short target regions [17,43]. The mitochondrial COI sequences of marine species belonging to various invertebrate taxonomic classes shared habitat with B. californica: Demospongiae of the phylum Porifera; Hydrozoa, Scyphozoa, and Anthozoa of the phylum Cnidaria, Ascidiacea of the phylum Chordata; Hexanauplia of the phylum Arthropoda; Crinoidea, Asteroidea, Ophiuroidea, Echinoidea, and Holothuroidea of the phylum Echinodermata; and Gymnolaemata of phylum Bryozoa were obtained from GenBank; accession numbers are presented in Table A1. These sequences were aligned using Clustal X [40], and analyzed to determine the regions conserved for B. californica, but sufficiently variable in related species. We selected the B. californica-specific nucleotide regions for primer and probe binding sites. Also, to increase the specificity for detecting target species, we designed taqman probe hybridized with only B. californica (Table 2).
To determine primer specificity, we used genomic DNA of species belonging to the classes of Hydrozoa, Ascidiacea, Crinoidea, Asteroidea, Ophiuroidea, Echinoidea, Holothuroidea, and Gymnolaemata ( Table 3). The B. californica-specific region was amplified using the same PCR mixture and thermal cycling conditions used for the amplification of COI, except for the primer (species-specific primers developed in this study) and annealing temperature (60 • C). The PCR products were separated by electrophoresis on a 1.5% agarose gel stained with ethidium bromide (Bioneer, Daejeon, Korea), and 2 µL of 100 bp DNA ladder (Elpis Biotech, Daejeon, Korea) was loaded to show the size of the proteins.
The B. californica-specific region was amplified by qPCR using a buffer solution containing 10 µL of qPCR BIO Probe Mix Hi-ROX (PCR Biosystems, London, UK), 5 µL of template eDNA, 1 µL of TaqMan probe at 5 pmol (Metabion, Martinsried, Germany), and 1 µL of each primer at 10 pmol (Cosmogenetech, Seoul, Korea) and 2 µL of DEPC-treated water. The cycling protocol, with optimal temperatures, was 95 • C for 3 min, followed by 45 cycles at 95 • C for 5 s (denaturation) and 60 • C for 30 s (annealing/extension), using an Applied Biosystems thermal cycler (Applied Biosystems, Waltham, CA, USA).    Table 1). The seawater was first filtered through a 300 µm nylon mesh for removing large clogs and then through a 3.0 µm nitrocellulose membrane (Merck Millipore, Darmstadt, Germany) [26,27,44,45]. The membrane filters were stored in a sample tube, and the sample tubes were stored in an ice box containing dry ice (−70 • C), moved to the laboratory as soon as possible, and then processed to isolate eDNA. eDNA was isolated from the membrane using a genomic DNA extraction kit (Bioneer, Daejeon, Korea) in accordance with the manufacturer's instructions, but with a modifications (400 µL ddH2O was used for elution instead of the elution buffer in the kit). The quality and concentration of eDNA were determined using a Nanodrop ND-1000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA), and all eDNA samples were stored at −80 • C until use.

eDNA Detection
The presence and biomass of B. californica DNA in the eDNA samples were analyzed by qPCR with species-specific primers and probe that were developed in this study. The probe was designed as a nucleotide sequence with high GC content and a high probability of binding (showing the low intraspecific variation or the conserved region in B. californica) in the sequence amplified with species-specific primers. qPCR was conducted using same qPCR mixture and thermal conditions used to confirm the sensitivity of the molecular makers with qPCR.
A B. californica-specific DNA fragment was used for DNA cloning (GNC Bio, Daejeon, Korea). DNA for cloning was subjected to 10-fold serial dilutions and a standard curve was plotted for quantification. We estimated the quantity of B. californica DNA in 8 L of seawater using a formula obtained from the y-intercept and slope, from which the number of detected DNA copies could be calculated. The number of copies of B. californica DNA was estimated by substituting the Ct value in the formula to estimate species biomass. PCR efficiency was calculated by substituting the slope value in the following formula: PCR efficiency (%) = ((10(−1/slope)) − 1) × 100% [46]. Also, we examined the limit of detection (LOD) and limit of quantification (LOQ) for the B. californica eDNA analysis as the lowest values of the linear range covered by the standard curve we established. A-cut-off values were calculated according to the method used by Kim et al. (2020) [26]. The reliable values in qPCR analysis were used as data and presented in figure.

Gene Selection for Designing Species-Specific Markers
A 658 bp region of B. californica COI was successfully amplified and sequenced with universal primers, and a 507 bp region of eight COI mt-DNA sequences was obtained from GenBank. There was no intraspecific variation within B. californica; the interspecific variation between the other four Bugulina species was 18.3-25.7% (Table A2). In the phylogenetic tree, B. californica from Korea formed a clade with the voucher specimen of B. californica and was clearly distinguished from the clade of related species (Figure 2).

Species-Specific Markers Design and Validation
The specific primer pair was designed based on the barcode region of mitochondrial COI of B. californica and 138 related species (Tables 2 and A1), namely BuCa_SF (forward), BuCa_SR (reverse), and BuCa probe (taqman probe). Designed species-specific primers within the amplified COI mtDNA region can hybridize only with the B. californica sequence at an annealing temperature of 60 • C, because of the differences in all the nucleotide sequences of the compared species, and the probe supported the detection of only the target species in qPCR, designed based on sequences amplified with species-specific primers. Template DNA from each species, including the various taxonomic species, were to confirm the specificity of the species-specific primer set, and a single, clear 122 bp band for each of three B. californica DNA samples was obtained on the agarose gel electrophoresis using the primer pair (Table 3, Figure 3). Subsequently, primer pair specificity was verified by PCR amplification for the different species, including native and non-native species inhabiting Korean coasts (Table 3). During specificity analysis of the primer pair using genomic DNA isolated from the other species, no amplification was observed ( Figure 3).
For develop more accurate and efficient marine bryozoan species monitoring method, we designed the specific molecular markers for B. californica. We produced the standard curve for our molecular markers. The number of DNA copies was estimated by a standard curve using 10-fold plasmid dilution ( Figure 4) and calculated using the formula obtained from the standard curve. The PCR efficiency calculated using the slope (−3.4257) was 96% (r2 = 0.9974).  Among the 12 sites investigated, B. californica colonies were detected in seven sites (1, 2, 6, 7, 8, 9, and 11); all of the survey sites showed amplification of B. californica DNA in all seawater sample replicates, whereas B. californica colonies were not detected at five sites (3, 4, 5, 10, and 12) using traditional survey methods (Figure 1). The Ct values for all qPCR reactions ranged from 23 to 44 and the negative control was confirmed in each experimental plate. We calculated the number of DNA copies using the y-intercept and slope after substitution of the Ct value into the formula (Figure 4).

Calculation of B. californica DNA Copies
The concentration of B. californica DNA in seawater samples between May 2017 and June 2018 is shown in Figure 5. B. californica DNA was detected in all survey sites, except one site (12); however, Ct values were not measured for several months in site 12. In the investigational period (August-November 2017, and March-May 2018) for which DNA was detected in all of the investigated sites, the amount of DNA was lowest in site 12. We calculated the average of the number of DNA copies per month during the whole investigation period, and the range was 2274-598,117 copies/L. The months were divided in three groups, depending on the number of DNA copies (copies/L): ≤10,000 (May, June, September, October, 2017, and March, 2018); 10,000-300,000 (July, November, 2017, and January, February, June, 2018); 300,000-600,000 (August, December, 2017 and April, May 2018). In particular, in August 2017, December 2017, and April 2018, more than 500,000 copies/L of DNA were detected, which was the highest number. There was no significant difference in the total amount of DNA detected between the sea areas, and the peak was slightly different by months at the survey sites; however, seasonally, the peaks were in late summer and early autumn (August-September), early winter (November-December), and spring (March-May). The amount of the detected DNA was divided into seasons and the averages were calculated: 211,787 copies/L in spring (March-May); 181,231 copies/L in summer (June-August); 8705 copies/L in autumn (September-November); and 240,025 copies/L in winter (December-February).

Discussion
The molecular markers for the detection of the invasive bryozoan, B. californica, designed in this study were sensitive and accurate. The specificity of molecular markers is important for the detection of species as there may be DNA of other species in mixed samples, such as environmental samples. When developing species-specific markers, it is important to identify the target species clearly and to find a suitable region of the gene for the preparation of markers. To detect only the target species in the environmental samples, the molecular markers used should be specific and not universal. Therefore, the DNA barcode region, with a large interspecific variation and a small intraspecific variation, is appropriate for the base sequence of marker design. In this study, we selected the universal animal DNA barcoding region (658 bp of COI) to design a species-specific marker for B. californica. Therefore, we verified species identification by performing molecular phylogenetic analysis based on COI sequences and confirmed that samples we collected in Korea formed a monoclade with the B. californica voucher specimen. In addition, it was found that non-intraspecific variation and large interspecific variation (18.3-25.7%) in bryozoan COI sequences indicated an appropriate genetic site for the design of species-specific markers ( Table A2).
The traditional method for species identification is based on morphological characteristics that may be affected by environmental factors, such as habitat and food, and relies on the traits of adult specimen; thus, skilled professional training and experience is required [47]. Furthermore, experts cannot often identify immature stages, juveniles, or broken specimens [13,48], and culturing the immature stages of organisms to adults for identification may be a lengthy process. These problems require a rapid and accurate new approach, such as species-specific markers that are designed based on specific nucleotides of target species. The specific genetic markers can hybridize with particular sequences only present in the target species, preventing misidentification and allowing the identification of all stages of the life cycle, even the larval stage [6].
We designed the B. californica-specific primers to amplify a 122 bp region of the mitochondrial COI sequence for eDNA analysis (Figure 3), and it is important to increase the possibility of target DNA detection in eDNA samples, which are degraded into short fragments by the influence of environmental factors, such as water, UV radiation, enzymes, and the activity of bacteria and fungi [49,50]. Several studies reported that eDNA is usually present as short DNA fragments (300-400 bp), and that a lake water temperature of 18 • C will maintain a 400 bp fragment of eDNA for approximately 1 week [51]. Therefore, Hajibabaei et al. (2006) [9] suggested that when detecting target DNA from degraded DNA samples, such as eDNA, designing species-specific markers that amplify only a short length are beneficial for preventing non-detection. In our study, the B. californica-specific primers and probe we developed were able to successfully detect the B. californica DNA in eDNA samples isolated from sea water (Figure 1).
We examined the distribution and abundance of B. californica in Korea using eDNA and compared the efficiency of eDNA analysis with the field survey. In Figure 1, more distribution sites were detected from qPCR analysis based on eDNA than field investigation that relied on visual evidence. B. californica DNA was detected at sites (3,4,5,10,12) where no evidence of B. californica was found in the field investigation. This suggested that B. californica is an aquatic organism that is difficult to detect by traditional methods, and showed that eDNA methods were effective for determination of the distribution of aquatic organisms [5,7] and were applicable to invasive bryozoan species in the marine environment. Therefore, our results may be due to the limitations of field survey methodologies that target aquatic organisms. However, in this study, we examined unpredictable sites using traditional methods to detect target species specific nucleotides from eDNA that was isolated from the imprints of organisms remaining in the environmental sources, without directly detecting the specimen; our data also suggested that eDNA methods can support field data and are applicable to marine environments.
In addition, we monitored monthly changes in the number of B. californica DNA copies detected in eDNA samples from May 2017 to June 2018 ( Figure 5). In our results, the amount of DNA was slightly different at each survey site, but peaked in late summer, early winter, and spring. Maturo (1959) [52] reported that B. californica present at Beaufort in North America reproduced from April to early December, and that colonies were most abundant during early September and December of 1954 and late April 1955. This period appeared to be similar to the peak of amount of DNA copies showed in our results ( Figure 5). Both laboratory and field studies have shown that an increase in abundance or the density of target species can lead to an increase in either eDNA concentration [16,18,45,53,54] or eDNA detectability [55]. Therefore, it can be suggested that eDNA derived from adults, rather than from larvae, appeared to have a greater effect on the results in this study. On the assumption that the adults are continuously present, if larvae, which are continuously released from late spring to early winter, had a larger effect on eDNA, larger amounts of DNA should have been detected during the release period of the larvae (larvae + adults) than in other periods (only adults). However, our results did not yield the expected results. Rather, the amount of DNA showed a difference when seasonally divided, and was highest in winter. It is expected to be the result of the seasonal difference in the water temperature affecting the maintenance of the eDNA in the environmental sample. In aquatic systems, environmental conditions influence eDNA persistence [49]. According to Strickler et al. (2015) [50], the colder the water, the greater the protection from solar radiation, and more alkaline environments are likely to hold detectable amounts of eDNA for longer than those that are warmer, sunnier, neutral, or acidic; in addition, warm water caused greater degradation of DNA by contributing to favorable environments for microbial activity. Thus, it can be supposed that the highest amount of B. californica DNA in winter was most likely due to greater preservation of eDNA as the water temperature is low. Nevertheless, the amount of DNA was lowest in autumn not summer. In summer, high temperature of seawater can reduce detectable eDNA templates by activating enzymes and microbial activities that lead to elevating the eDNA degradation rate. However, the B. californica has shown great abundance of colonies in summer, likely because B. californica has been reported to be more heat-resistant than other invasive bryozoan species and is less influenced by water temperature [56]. For this reason, the amount of DNA detected may not be the lowest in summer. Also, eDNA detection rate is related to the eDNA shedding rate of the target species. According to previous study, various abiotic factors-such as feeding activity, presence of co-habitat species, and temperature variation-can affect to eDNA shedding rate and degradation rate [57]. They suggested the presence of filter feeding organism as co-habitat species can cause the loss of the release DNA and affect to target species eDNA detection. Furthermore, eDNA degradation is faster in freshwater than seawater, and these showed eDNA sustained rate can increased by salinity variation with sea area [58]. However, all of these concepts were not been perfectly applied in various taxonomy and environment where many of co-habitat species exist. Each of species have different eDNA shedding rates, and there are so many types of ecosystem which have various abiotic factors and species composition. Therefore, it is expected the understating eDNA experimental study with many situation including various variables have a strong impact on future eDNA analysis.
During the investigation period in which DNA was detected in all surveyed sites, B. californica DNA was not detected in May−July 2017, October 2017, and June 2018 at site 12. It may be that the amount of target species DNA in environmental samples is below the threshold that can be detected using eDNA methods [15], and the causes of the absence of the B. californica colonies in the site 12 are thought to be invisible aquatic species and low population. The enrichment of eDNA samples may help to exceeding the threshold of detectable DNA copies, but [20] reported that altering the DNA concentration may affect PCR inhibitors and lead to false results of eDNA detection. However, if the concentration of DNA samples was low, increasing the number of PCR cycles may improve detection [59] and should be accompanied by sequencing analysis to confirm the amount and complexity of nonspecific background products [60,61].
These results suggested that the method based on the detection of DNA in environmental sources, combined with species-specific molecular markers, could be used to monitor the distribution and aspects of biomass changes as a complementary technique to traditional field survey methods.

Conclusions
In conclusion, application of environmental DNA (eDNA) method for the invasive species has invisible life stage can help the detection and monitoring. This method is efficient for the invasive bryozoan species, B. californica, in marine environments, as the species is difficult to detect using traditional methods, present at low densities, and at cryptic life stages, such as larva. Monitoring of the monthly changes in DNA concentration in eDNA can also be used to track increases or decreases in the abundance of the species and how various environmental factors affect the abundance of the target species, and offers advantages for managing the distribution and predicting diffusion. When the traditional survey was only method for monitor invasive bryozoan, B. californica, the colonies were not found in several investigation sites. On the other hand, when the molecular survey using eDNA method was applied, we can detect more B. californica distribution sites than colonies found sites through B. californica DNA detection in environmental samples. Additionally, this method make possible to examine the biomass of B. californica by calculation the DNA copies in eDNA detected sites. Our newly designed species-specific markers used in this study were shorter and B. californica specific, and it has advantage for the degraded and various species DNA mixed environmental DNA samples. Therefore, the eDNA method in this study can contribute for the supplement the limit of traditional survey, especially in case of marine invasive bryozoan species.