First Record Sparus aurata Larvae in Teluk Penyu Beach Cilacap, Indonesia, revealed by DNA Barcoding

Morphological identification often fails to reveal the taxonomic status of fish larvae. DNA barcoding using particular DNA segment is frequently successful in solving the faulty and might reveal overlooked species, including species outside their recognized geographic ranges, such as Eastern Atlantic seabream Sparus aurata. This study aims to assess fish larvae diversity in Teluk Penyu Beach, Cilacap Central Java, Indonesia, through the cytochrome c oxidase 1 barcoding. Fish larvae were collected using a larva net with a diameter of 60 cm and height of 150 cm and horizontal towing during the field trips in March and April 2021. Larvae morphotypes were determined based on their general morphological performance observed under a magnification lens according the available references. Eighteen larvae morphotypes were successfully barcoded, and 5% genetic distance was used as a species border. Fourteen fish larvae species were revealed, with intraspecific genetic distances between 0.00% and 4.12%, while intrageneric genetic distances ranged from 5.50% to 19.29%. An interesting finding was that one larva morphotype was barcoded as S. aurata with high genetic identity (99.19% to 99.68%) and low genetic distance (0.32% to 0.81%). The discovery provides the first new data on S. aurata in Teluk Penyu Beach Cilacap, Central Java, outside its recognized geographic range. Overall, this study provides information about fish larvae in Teluk Penyu Beach, essential for estimating the number and species of fish recruited in the southern Cilacap, which is vital for fisheries management. Nevertheless, new research requires further clarification of S. aurata in Teluk Penyu Beach.


Introduction
Sparus aurata Linnaeus 1758, known as Gilthead seabream, is well-known as a Sparidae member.It is an economically important marine fish in Central Europe and the Mediterranean Sea (FAO, 2006;Rossi et al., 2006;Arabaci et al., 2010).S. aurata had been proposed as a synonym of Pagrus auratus (Silver seabream) (Froese and Pauly, 2023; World Register of Marine Species).At present, the synonymy between S. aurata and P. auratus is rejected and referred to as misapplied names (Paulin, 1990).The rejection of synonymy between S. aurata and P. auratus) was supported by previous molecular studies, showing that both morphospecies are two different genetic species (Somarakis et al., 2013;Landi et al., 2014;Bernard-Capelle et al., 2015;Fotedar et al., 2019).
Six species of Sparidae have been described in Indonesia (Froese and Pauly, 2023).The number of species can be different from one to the other references.One logical argument is that a lot of synonyms are available.For example, Acanthopagrus pacificus (Frose and Pauly, 2023; World Register of Marine Species) has Acanthopagrus berda as one synonym, but a study uses Acanthopagrus berda (Limmon et al., 2019).None of the six sparid species in Indonesia belong to the genus Sparus (Froese and Pauly, 2023; World Register of Marine Species).Wibowo et al. (2018) reported Sparus larvae from the north coast of Jakarta, Indonesia, revealed from DNA barcoding.However, Wibowo et al. (2018) did not show to which species those larvae belonged to.Using a similar method and sample types (fish larvae), Azmir et al. (2017) reported Sparus aurata in the mangrove ecosystem in Malaysia.However, no study reported the presence of S. aurata in the southern Cilacap Coastal areas, especially in Teluk Penyu Beach.
Due to limited characteristics for precisely identifying fish larvae at the species level (Ko et al., 2013), exploring fish larvae diversity is still an exciting challenge for biologists and ichthyologists in Indonesia.This challenge has been successfully overcome by the development of molecular identification (Hubert et al., 2010).Numerous studies have proved that molecular identification, especially cytochrome c oxidase 1 (COI) barcoding, is reliable for fish larvae identification (Wibowo et al., 2015;Wibowo et al., 2016) and revealing undescribed and overlooked species worldwide (Baldwin and Johnson, 2014;Bhagawati et al., 2020;Brinocoli et al., 2020;Nuryanto et al., 2020).Besides revealing the described fish species from a particular site, DNA barcoding on fish identification might reveal undescribed and overlooked fish species, which was not revealed by morphological identification (Wibowo et al., 2018;Winarni et al., 2023).
Coastal areas of southern Java Island are directly connected to the Indian Ocean.This geographical position allows southern Java to receive fish larvae from the Indian Ocean through ocean currents.However, more studies about fish larvae diversity in the southern coastal areas of Java need to be done.Only one study reported larvae diversity from southern coats of Java through DNA barcoding (Nuryanto et al., 2017).In their research, Nuryanto et al. (2017) collected larvae from the eastern areas of the Segara Anakan estuary in the Cilacap Regency, Central Java.The sampling locations were gone further inside the estuary so that it could not describe the diversity of fish larvae on the south coast of Cilacap.Therefore, data on fish larvae from the south coast of Cilacap still needs to be included.This study aims to assess fish larvae diversity in Teluk Penyu Beach, Cilacap Central Java, Indonesia, through the cytochrome c oxidase 1 barcoding.The data have important implications for further study to ensure that S. aurata lives in Indonesia's marine ecosystem and the management of the marine fisheries, especially on the southern coast of Cilacap, Central Java, Indonesia.

Materials and Methods
Fish larvae were collected at Teluk Penyu Beach Cilacap Regency, Central Java, Indonesia.This study retrieved a map to show sampling sites from the online version of Google Earth in September 2021.The sampling tracks in Teluk Penyu Beach are illustrated in Figure 1.

Field trips
Field trips for fish larvae sampling was conducted in March and April 2021.Larvae were collected using a larva net with a mouth diameter of 60 cm and trapezium height of 150 cm.At the rear end of the trapezium is attached the sample bottle.We conducted sampling in the morning and evening (Nuryanto et al., 2017) to ensure that diurnal and nocturnal larvae were collected.Larvae were sampled in seven tracts at the study site (Figure 1).Four-time horizontal towing was conducted in each sampling tract; therefore, the total efforts were 28 towing.Towing was conducted by tightening larva nets in the rear part of a boat, and the driving speeds were approximately 3 knots.The time duration for each towing was 10 minutes.After ten minutes of driving, the larva net was filled up, and the collected materials inside the sample bottle were poured into a flour sieve (Nuryanto et al., 2017;Kusbiyanto et al., 2020).Collected materials were soaked in 96% technical ethanol and washed in 70% alcohol to quickly sort fish larvae and separate them from larvae of other organisms and garbage.Fish larvae were put in new sample bottles filled with fresh 70% ethanol for temporary preservation.Samples from each towing were placed in different bottles and labeled.

Larvae sorting and permanent preservation
Upon arrival in the laboratory, fish larvae samples were then washed and sorted carefully to avoid the presence of larvae of other organisms and debris still mixed.Clean and good conditions larvae were then preserved permanently in new 70% ethanol for examination.

Morphotype identification
Instead of having a transparent body, fish larvae samples had white body color after ethanol preservation.That larvae condition made it easier to differentiate among morphotypes under a magnification lens.Morphotype identification followed Leis dan Carson-Ewart (2000) and Leis (2015).Examples of larvae morphotypes are presented in Figure 2.

Larvae barcoding process
One specimen from each morphotype was shipped to a company for DNA barcoding.Genome was extracted using DNA Miniprep Kit (Zymo Research, D6016) following the manufacturer's protocol.The Polymerase Chain Reaction (PCR) amplification of the COI gene was performed using the MyTaq HS Red Mix (Bioline, BIO-25047) and universal primer pair FishF2-5' TCGACTAATCATAAAGATATCGGCAC3' and FishR2-5' ACTTCAGGGTGACCGAAGAATCAGAA3' (Ward et al., 2005).The amplification settings were started with an initial denaturation at 96° for 3 minutes.The process was continued with denaturation at 94° for 10 seconds, annealing at 52° for 30 seconds, and extension at 72° for 45 seconds with total cycles were 35 times.The volume of each chemical component for the final volume 25 µl PCR mixtures was KOD FX Neo 1 µl, 2X PCR Buffer KOD FX Neo 12.5 µl, 2mM dNTPs 1 µl, 10pmol/µl of each primer was 1 µl, template DAN 1 µl, and ddH2O 6 µl.Total genomic DNA isolation and COI marker amplification were conducted at Genetika Laboratory (PT.Genetika Science Indonesia, Jakarta).The sequencing of the COI gene was used as the bi-directional sequencing technique at 1 st BASE Asia in Kuala Lumpur, Malaysia.

COI sequences editing and trimming
Either forward or reverse sequences of the COI gene from each morphotype were selected for further analysis to have homogenous data.This study used one sequence direction because some morphotypes had a mesh-up chromatogram.The fragments of the COI gene were edited and trimmed in BioEdit 7.0 (Hall, 2005).A functional portion of the COI was tested by translating the obtained sequences into amino acid sequences using an online version of ORF finder (https://www.ncbi.nlm.nih.gov/orffinder) and through basic local alignment search tool (BLAST) alignment to check whether a gap was present or not along with the sequences.

Barcoding Analysis
Multiple pair sequence alignment was conducted using ClustalW in BioEdit7.0(Hall, 2005).The present study utilized sequence identity, genetic distances, and monophyly data to delineate taxonomic status.A homology test was performed through BLAST comparison to the closest relative in the GenBank.An identity value of 96% was chosen as the identity threshold for species level (Mat Jaafar et al., 2012Pereira et al., 2013), considering the geographic locality of the sample and reference species (Higashi et al., 2011).Sequence similarity Figure 2. Some of the fish larvae collected during the samplings was double-checked with the nearest references in BOLDSystem (Ratnasingham and Hebert, 2007).In the case of high genetic identity and multiple closest reference species were observed in both databases.This study also considered the reference species' expected value and geographic range during taxa delineation.
The Kimura 2-parameter (K2P) genetic distances were calculated between samples and the nearest reference species.In line with 95% sequence homology, a 5% genetic divergence was defined as a genetic gap within species because each fish species has variable genetic distance values (Pereira et al., 2013;Lim et all, 2016).
This study used Lates calcarifer as an outgroup species to obtain reliable polarity of the tree.The support to sequence homology and genetic distances came from the monophyly of samples and nearest relative in GenBank.The monophyly samples were obtained from the K2P neighbor-joining (NJ) tree and reconstructed in MEGA XI (Tamura et al., 2021).Tree reconstruction was repeated for s 1000 bootstraps replication.
The study also calculated the K2P genetic distance between TP06 and two species of Sparidae (S. aurata and P. auratus).Based on previous literature, we chose those two synonymized species (Fricke et al., 2021).The phylogenetic between TP06 and reference species was also reconstructed using the neighbor-joining algorithm and K2P substitution model in MEGAXI (Tamura et al., 2021).Lates calcarifer was used as the outgroup species.The bootstraps value obtained from 1000 bootstraps replications also supported the branching pattern.

Result and Discussion
A total of 632 larvae of fish species were obtained during the sampling at the study site.Based on the body's general shape, 19 larvae morphotypes were observed.Homology test resulted in genetic identity ranging from 83.42% to 100% to the nearest five top hits conspecific reference in GenBank with expected value going between 0.00 for most reference species and 4 -176 for TP02 to B. dussumieri (Table 1).After several tries, TP10 was still not successfully barcoded.Therefore, only 18 of 19 morphotypes were gained and utilized in further analysis.
Previous studies proved that the natural population shows wide genetic divergences (Mat Jaafar et al., 2012;Pereira et al., 2013;Lim et al., 2016;Nuryanto et al., 2018).The genetic value ranges from 3% to 5%, referred to as a moderate value, and should consider additional information (Jeffrey et al., 2011), such as geographic locality (Higashi et al., 2011;Candek and Kuntner, 2015).Candek and Kuntner (2015) stated that 5% genetic divergence is a moderate value.This study considers geographic locality between our samples and references species, especially for genetic divergent above 3% (genetic identity between 96.01% and 96.37%).For example, reference species of morphotype TP015 (C.bilineatus) were collected from Taiwan (Chang et al. 2017) and the South China Sea, thousands of kilometers away from Teluk Penyu Beach Cilacap, Central Java, Indonesia.It is reasonable that this study used 5% genetic divergence as the threshold value for species delineation.Some other morphotypes could only be assigned at the genus level because their homology value was below 96% (between 83.42 and 94.74%) and delineated into five species.Those morphotypes were TP03 (Parioglossus sp.), TP08 (Stolephorus sp.), TP09 (Amoya sp.), TP13 (Stolephorus sp.), TP17 (Eugnathobius sp.), TP18 (Redigobius sp.), and TP19 (Amoya sp.).The assignment of those seven morphotypes into genus levels was convincing and reliable.The placement of those morphotypes based on sequence divergences ranged between 5.25% and 16.58% to their reference species.Zhang and Hanner (2011) reported that sequence divergence among congeners was 17.16%.This study obtained genetic divergence among congeners below the value reported by Zhang and Hanner (2011).Intraspecific K2P genetic distance observed in this study ranged from 0.00% to 4.12%.The maximum K2P genetic distance was below 4.8% reported for Genus Hemiramphodon (Lim et al., 2016) and Carangidae (Mat Jaafar et al., 2012).Even the maximum intraspecies K2P genetic distance obtained in this study was below 8.5% in Neotropic fish (Pereira et al., 2013).Therefore, the delineation of TP15 as the larvae C. bilineatus was reasonable.This delineation was supported by the geographic distribution of C. bilineatus, which includes the Indo-Pacific region (Menon, 1977;Froese and Pauly, 2023).Intrageneric K2P genetic distances ranged from 5.50% to 19.29%, with a mean value of 14.47%.The maximum value was below the maximum value reported by Pereira et al. (2013), which was 24.8% for Neotropic fish, and Ahmed et al. (2019), who wrote maximum intrageneric K2P genetic distance was 27.03% for indigenous fish in Bangladesh.The mean value was also below 18.25% for ophichthid fishes (Xing et al., 2020) and engraulids fish (Afrand et al., 2018).
The morphotypes TP07 and TP20 were barcoded as single genetic species (Ambassis sp.).The present study delineated morphotypes TP04 and TP16 as N. soldado.We found a similar result for morphotypes TP08 and TP13, identified as Stolephorus sp.The finding was logical because fish larvae from single species might have different morphology at different stages (Ko et al., 2013).In contrast, fish larvae of different species might share similar morphology (Victor et al., 2009).A previous study reported a similar result on fish larvae from eastern areas of Segara Anakan estuary, Cilacap (Nuryanto et al., 2017).
Five top hits specimens in GenBank for TP06 consisted of four S. aurata and one Omobranchus ferox with high identity, expected' value of 0.00, and low genetic distances (Table 1).S. aurata and O. ferox are distantly related species because they belong to different orders.S. aurta is from Eupercaria, while O. ferox belongs to Blenniiformes (Froese and Pauly, 2023).Further confirmation was performed to ensure the taxonomic status of TP06 by doing a similarity check to data in BOLD (Ratnasingham and Hebert, 2007).The result showed that over 10 top hit specimens for TP06 in BOLD were S. aurata, and none of them was O. ferox.Even after the top hits search was increased to 100 samples, no O.ferox was observed.This study conducted a further analysis by comparing it to S. aurata and Pagrus auratus to be more convinced that the Morphotype TP06 is Sparus aurata larvae.Both species were compared because both scientific names had been synonymized (Froese and Pauly, 2023;Fricke et al., 2021).The average K2P genetic distance between morphotype TP06 and two sparid species was 0.47% to S. aurata and 21.80% to P. auratus, respectively (Table 2).
An additional phylogenetic relationship between morphotype TP06 and two sparid species was conducted (Figure 3).It is seen in Figure 3 that morphotypes TP06, S. aurata, and P. auratus formed a monophyletic clade (A: Sparidae).That clade separated from the outgroup specimen (Lates calcarifer GU674017, Latidae).Sparidae clade was divided into subclades B and C. Subclade B was formed by morphotype TP06 and S. aurata specimens.This subclade was supported by the highest bootstrap value (100/100) for ML and MP algorithms with branch lengths less than 0.020 (Figure 3).Subclade C consisted of all five P. auratus retrieved from GenBank.Based on genetic identity, distance, and evolutionary relationships, this study finally defined TP06 as genetically identified as Sparus aurata.
The assignment TP14 into a species level based on genetic identity, expect value, genetic distance, and monophyly.Morphotype TP14 has two nearest reference species in GenBank: Favonigobius gymnauchen and Oxyurichthys microlepis.A similar result was obtained when a similarity test was conducted on the data in BOLD (Ratnasingham and Hebert, 2007).The result still needed to be clarified.Morphotype TP14 was further analyzed by comparing it to five species in GenBank that showed the highest identity (Table 3).All five species retrieved from GenBank belong to Gobiidae.The K2P genetic distance between morphotype TP14 to five gobiid species ranges between 0.65 and 17.57 (Table 3).The K2P genetic distance comparison to five gobiid species indicates that TP14 is genetically closer to F. gymnauchen than O. microlepis (Table 3).To further improve the reliability of TP14 placement to F. gymnauchen, we analyzed the phylogenetic relationship between TP14 and five gobiid species (Table 3).The phylogenetic tree was reconstructed using maximum likelihood (ML) and maximum parsimony (MP) algorithms (Figure 4).We could infer from Figure 4 that morphotype TP14 formed a monophyletic subclade with F. gymnauchen (A).Morphotype TP14 has a closer evolutionarily relationship to F. gymnauchen than O. microlepis.Therefore, we could assign that TP14 belongs to F. gymnauchen.
Sparus aurata is a seabream species distributed in the Eastern Atlantic and Mediterranean Sea (Rossi et al., 2006;FAO, 2006;Arabaci et al., 2010;Bodur, 2018).No study recorded the existence of Sparus aurata from Indonesia (Froese and Pauly, 2023).Nevertheless, official report of Fish Quarantine and Inspection Agency (FQIA), Ministry of Marine Affairs and Fisheries, Republic of Indonesia reported that S. aurata species already exist in Indonesian waters.However, BKPIM did not explain the exact location where S. aurata was found in Indonesian waters (http://www.bkipm.kkp.go.id/ bkipmnew_rubah/ias).
S. aurata has been synonymized with P. auratus.However, the result of this study did not support the synonymy between S. aurata and P. auratus.All the data in this study strengthen that S. aurata is different from P. auratus, described by Paulin (1990).Dahruddin et al. (2016) reported one species of Sparidae from Java but did not explain its species.Therefore, this study provides the first data about S. aurata on the southern coast of Cilacap, Central Java, Indonesia.This information at least has two meanings.First, the geographic range of S. aurata is possibly much broader than previously suggested (FAO, 2006;Rossi et al., 2006;Arabaci et al., 2010).Second, S. aurata larvae have been transported from the Mediterranean Sea to Teluk Penyu Beach Cilacap through water ballast from marine shipping.
Transport through shipping is possible because S. aurata has a long planktonic larval duration of approximately 43 and 50 days (Froese and Pauly, 2023).
The existence of S. aurata outside of their natural range has been reported by Azmir et al. (2017) in the mangrove ecosystem in Malaysia.Other species of Sparus, i.e., Sparus sp., have been reported from coastal areas near Jakarta (Wibowo et al., 2018).The present study and those previous studies (Azmir et al., 2017;Wibowo et al., 2018) obtained data from larvae barcoding.However, no data available about S. aurata in Indonesia got from the study on the adult individual.The existence of fish species outside their geographic ranges was reported by Muchlisin et al. (2017), who found Anguilla bengalensis in Aceh waters, Indonesia.Anguilla bengalensis was a well-known endemic in Bengala Bay, India.
Further study on S. aurata in Indonesia is highly needed, emphasizing the adult individual.Comprehensive data about S. aurata in Indonesia, covering all life stages, have important implications for the control of S. aurata in Indonesia marine water since it has been listed as invasive species (http://www.bkipm.kkp.go.id/bkipmnew_rubah/ias).
A total of 14 taxa were obtained during this study.That number was higher than that reported by Nuryanto et al. (2017), who collected fish larvae from the eastern areas of the Segara Anakan estuary.On the one hand, the present study reported similar fish larvae with the survey by Nuryanto et al. (2017).Both studies found F. gymnauchen and Stolephorus sp.In

Conclusion
Fish larvae from Teluk Penyu Beach Cilacap were successfully barcoded and placed into 14 taxa.Sparus aurata was among the fish larvae found during the study.The existence of S. aurata in Teluk Penyu Beach Cilacap suggested that the geographic distribution of S. aurata is much broader than previously reported only in Eastern Atlantic and the Mediterranean Sea.Current research findings are vital for sustainable fisheries management in the Cilacap Region in Central Java, Indonesia.This study only sampled specimens from Teluk Penyu Beach; only limited morphotypes were sequenced.Representatives from other locations and more sequenced samples are needed to provide statistically reliable data for fisheries and conservation.

Figure 1 .
Figure 1.Sampling site showing sampling tracts during larvae collection (red lines: sampling tracts)

Figure 3 .
Figure 3.The ML and MP trees show monophyly between TP06 and S.aurata with high bootstrap support (left: MP value; right: ML value).

Table 1 .
Genetic identity and K2P distance between each morphotype and their conspecific references

Table 1 .
Genetic identity and K2P distance between each morphotype and their conspecific references (continue)

Table 3 .
The K2P genetic distances between morphotype TP14 and five gobiid species of the remaining species were different.Therefore, the present studies increase the data on fish larvae diversity in Cilacap marine waters, Central Java, Indonesia.Moreover, despite the limited number of larvae obtained and the narrow sampling location, this study provides the first and new insight into the presence of S.aurata in Indonesian waters, particularly in Teluk Beach Cilacap, Central Java.These data are needed to estimate the number and species of fish recruited in the southern Cilacap, which is essential for fisheries management.