Phylogenetic analyses reveal bat communities in Northwestern Mexico harbour a high diversity of novel cryptic ectoparasite species

Background Parasites are integral parts of ecosystem function and important drivers of evolutionary processes. Characterising ectoparasite diversity is fundamental to studies of host-parasite interactions, evolution, and conservation, and also for understanding emerging disease threats for some vector borne pathogens. With more than 1,400 species, bats represent the second most speciose mammalian clade, but their ectoparasite fauna are poorly known for most species. Methods We sequenced mitochondrial Cytochrome Oxidase C subunit I and nuclear 18S ribosomal gene fragments, and used Bayesian phylogenetic analyses to characterise ectoparasite taxon identity and diversity for 17 species of parasitized bats sampled along the Baja California peninsula and in Northwestern Mexico. Results The sequence data revealed multiple novel lineages of bat bugs (Cimicidae), flies (Nycteribiidae and Streblidae) and ticks (Argasidae). Within families, the new linages showed more than 10% sequence divergence, which is consistent with separation at least at the species level. Both families of bat flies showed host specificity, particularly on Myotis species. We also identified new records for the Baja peninsula of one tick (Carios kelleyi), and of five Streblid bat flies. One Nycteribiid bat fly haplotype from Pallid bat (Antrozous pallidus) hosts was found throughout the peninsula, suggesting potential long distance co-dispersal with hosts. Different bat bug and tick communities were found in the north and south of the peninsula. Conclusions This study is the first systematic survey of bat ectoparasites in the Baja California peninsula, discovering highly genetically differentiated novel lineages compared to other parts of North America. For some ectoparasite species, haplotype distributions may reflect patterns of bat migration. This work is a first step in characterising ectoparasite diversity over the Baja California peninsula, and understanding how ecological and evolutionary interactions shape bat ectoparasite communities among host species in different parts of their ranges.

were identified by molecular assays (Najera-Cortazar, 2020). A list of sites and bat species 130 sampled is given in Table 1. Ectoparasites were collected manually from bats using forceps 131 and transferred to labelled Eppendorf tubes with 96% ethanol for storage. Ectoparasites 132 collected from the same bat but from different taxonomic families were stored in separate 133 tubes for each family, with corresponding labels and appropriate specimen source identifiers. 134 All samples were kept on ice during fieldwork, until their arrival to the laboratory where they 135 were stored at -20° C. Ectoparasites were photographed during fieldwork using a portable 136 Maozua 5MP 20-300X USB microscope, taking ventral and dorsal views. Thermal cycling parameters were as follows on a TECHNE thermocycler model TC-512: 159 initial denaturation step at 95°C for 5 min, followed by 40 cycles of denaturation at 94°C for 160 40 seconds, annealing at 53°C -56°C for 1 minute and extension at 72°C for 1 minute. Final 161 extension was performed at 72°C for 10 minutes. 162 For the 18S rDNA gene, approximately 800bp amplicons were amplified using the 163 primers 18S-1F (5' CTGGTTGATCCTGCCAGTAGT 3') and 18S-3R (5' 164 8 GCATCACAGACCTGTTATTGC 3') for flies (Whiting, 2002); and D-F (5' 167 GGCCCCGTAATTGGAATGAGTA 3') and C-R (5' CTGAGATCCAACTACGAGCTT 3') 168 for ticks (Mangold, Bargues, & Mas-Coma, 1997). The same reaction mix quantities were 169 used as for the COI gene, with MgCl2 at a final concentration of 1.5 mM for all primer sets. 170 Thermal cycling parameters were: initial denaturation at 95°C for 5 minutes, followed by 30 171 cycles of denaturation at 94°C for 40s, annealing at 53°C for 1 min and extension at 72°C for 172 1 minute. Final extension was performed at 72°C for 10 minutes. 173 PCR products were visualized by gel electrophoresis using 1% agarose with GelRed® 174 (Biotium Ltd) staining, and then sent for PCR product purification and Sanger sequencing to 175 Genewiz Inc. (Azenta Life Sciences, Leipzig, Germany), with each amplicon sequenced in 176 both forward and reverse directions. 177

Quality control and references sequences 178
Sequence quality for both the forward and reverse strands of each amplicon was 179 evaluated in BioEdit 7.2.5 (Hall, 2005). Trimmed forward and reverse sequences were 180 combined to generate a consensus sequence for each amplicon, and then analysed in BLAST 181 (The National Library of Medicine, 2018) to generate initial taxon identities and identify 182 reference sequences. Additionally, COI sequences were also analysed in the BOLD system 183 platform (Ratnasingham & Hebert, 2007   For COI, Cimex 1 is represented by a single specimen (EBCO155) obtained from a 274 Myotis yumanensis host at Mosqueda, and forms a sister lineage to Cimex 2, represented by 275 two haplotypes from two bugs parasitizing M. californicus hosts, which were also sampled 276 from northern sites. The 18S analysis also placed EB18S155 in a separate clade, but Cimex 2 277 is paraphyletic. Both markers placed Cimex latipennis, which is distributed from Canada to 278 the north western USA (Usinger, 1966), as the closest named molecular reference species to 279 these 2 lineages. Cimex 3 is represented by a monophyletic clade for COI, consisting of two 280 haplotypes, sampled from specimens collected from four Antrozous pallidus and one M. 281 californicus hosts, all distributed in the southern half of the peninsula (sites Rosarito, San 282 pilosellus and C. brevis formed a clade which appeared to be ancestral to Cimex 3, with 13% 285 divergence (Table 3). There is no molecular reference for C. brevis in the 18S analysis, and 286 Cimex 4 is the closest sister lineage. 287 The Cimex 4 COI lineage includes 17 haplotypes derived from 22 specimens, where 288 21 of the bugs were found parasitizing Parastrellus hesperus individuals, and one from A. incrassatus is reported as occurring on A. pallidus hosts (Usinger, 1966), but no reference 295 sequences are available. Therefore, one of the novel bug lineages may represent this species, 296 but further morphological assessments are required for confirmation. 297

Phylogenetic assessment of bat fly sequences 298
The best fit sequence evolution models were GTR+G for the COI gene set, and 299 T92+G+I for the 18S data. In total, 77 bat flies were sequenced, yielding 76 sequences for 300 COI and 73 for 18S amplicons respectively. Individual sequences for two specimens for COI, 301 and three for 18S did not pass sequence quality thresholds and were discarded. In the final 302 sequence set for the COI marker, there were 49 sequences from the Nycteribiidae family 303 (wingless bat flies) and 27 for the Streblidae family (winged bat flies), representing ten novel 304 lineages and six new species records for Mexico. For the 18S marker, 46 sequences were 305 generated for the Nycteribiidae family, and 26 sequences for the Streblidae family. For this 306 marker, nine novel lineages were observed (18S sequences corresponding to the Basilia 2a In the COI phylogenetic analysis, the 49 Nycteribiidae sequences formed five 309 lineages, all of which appeared to be novel with respect to GenBank references. Genetic 310 divergence among Baja lineages ranged from 2.9 to 14.5%, and up to 16% against reference 311 sequences (Appendix 5). For COI, Nycteribiid 1 and 2, formed sister clades with 4.4% 312 divergence, and 10% divergence from the closest reference haplotypes derived from Asian 313 Nycteribia, species. For 18S the closest references to Nycteribiid 1 and 2 were North 314 American species not available for COI, Basilia corynorhini and Basilia forcipata, 315 respectively. Nycteribiid 1 and 2 lineages were primarily associated with Myotis bat hosts, 316 but with one Nycteribiid 2 haplotype recovered from a Parastrellus hesperus host in the 18S 317 dataset (specimen EN18S514, Figure 4). 318 The three other Baja lineages formed clades associated with Basilia reference 319 sequences from species recorded in Madagascar, USA and Panama. Basilia 1 has 4.9%  and Artibeus hirsutus, respectively (Table 2). Most streblid lineages were parasitizing fruit-357 nectar feeding bats (Phyllostomidae) over the mid and northern peninsula, with the exception 358 of Trichobius 1 found on Mormoops megalophylla hosts, (family Mormoopidae), which feeds 359 on insects. Trichobius 1 was the only streblid fly lineage found in the south of Baja. The 360 other fly lineages were distributed in the southern continental sites (Tucson and Primavera 361 sites, Figure 1, Table 1). Outside of the closest sequence matches, none of the other species 362 reported as parasitizing hosts in Baja (Table 2) appear to have reference sequences in 363 GenBank. Genetic divergence within groups from the lineages in this study showed levels 364 between 0.0% and 1.1%, with the exception of A. phyllostomatis, which showed a 3.3% 365 divergence within its group (Appendix 6). 366

Phylogenetic assessment of Baja peninsula bat tick sequences 367
The best fit sequence evolution models were GTR+G+I for the COI gene set, and 368 K2+G+I for the 18S gene set. There were 45 ticks sequenced from the Argasidae family (soft 369 ticks), with 39 sequences generated for COI, and 44 for 18S, where six specimens failed to 370 amplify for COI (specimens ETCO157, ETCO306, ETCO338, ETCO480b, ETCO480c and 371 ETCO485), and one specimen failed to amplify for 18S (specimen ET18S457). One new 372 Baja species record and seven potential novel lineages were obtained ( Figure 6). The only 373 match with GenBank and BOLD COI records (95.95% and 97.07%, respectively) was the 374 soft tick Carios kelleyi, for 11 specimens (parasitizing A. pallidus hosts) with intra-clade 375 divergence of 1.94% to 2.47% ( Figure 6, Appendix 7). For COI data, Antricola 1 lineage had 376 10.6% divergence from Antricola marginatus, and 14.6% from A. mexicanus (Appendix 7). 377 Ornithodoros faccini had a divergence of around 11% for clades Carios kelleyi, Tick 1 to 5, 378 and O. yumatensis showed divergence of around 9.9% with respect to Tick 6 lineage. Ticks 379 of lineage Tick 6 were recovered from M. peninsularis from which there are no previous 380 records of bat ticks to our knowledge. 381 All sequences from this study generated the same lineages for both COI and 18S 382 genes ( Figure 6). Lineages C. kelleyi and Tick 1 to 5 formed a monophyletic group with 383 respect to the reference sequences with posterior support greater than 0.6 for both markers, Argasid ticks have previously been reported for the bat species in this study, primarily 486 Ornithodoros species (Table 2), but with the exception of Carios (Ornithodros) kelleyi, none 487 are close molecular matches for our sequences. In ticks, for COI, the threshold for between 488 genus divergence is considered to be above 10% (Hebert et al., 2003). The observed COI 489 divergence among lineages of this study (6.1% up to 19.3%), and to reference sequences 490 (9.9% to 22.8%), suggests that the Baja lineages could represent novel species and potentially 491 novel genera. 492 We found apparent host-specific and generalist lineages for the ticks reported here. 493 Lineages that appeared host-specific were restricted to single sites, while generalists were 494 found across multiple sampling locations. For example, Tick 3 parasitizing M. yumanensis 495 was found only in San Basilio, and Tick 6 was found only on M. peninsularis sampled in 496 Teste. Carios kelleyi was found to only parasitise A. pallidus in this study, and was found at three sites in the mid and north peninsula. However, C. kelleyi is known to parasitise multiple 498 species across North America and Cuba (Gill et al., 2004), and the reference COI sequence 499 used here derives from an Eptesicus fuscus bat sampled in New Jersey, eastern USA 500 (GenBank accession code: MT780277). Argasid ticks show a continuum of hosts-specificity 501 (Cumming, 1998;Esser et al., 2016), but tick stage-cycle must be considered, with immature 502 ticks being more generalist than adult conspecifics (Esser et al., 2016; Nava & Guglielmone, 503 2013). In this study, life stage was not assessed while collecting ticks, therefore it is possible 504 that there are gaps in host range regarding unidentified larvae that were not sequenced. 505 The Tick 5 lineage was found on multiple species in the mid and north peninsula, but 506 the only specimen recorded for the south peninsula was collected from Parastrellus hesperus 507 (specimen 338, Faro site, 18S only). Since P. hesperus is rare in southern Baja, the presence 508 of this tick indicates potential dispersal of P. hesperus from the north to the south peninsula. 509

Limits on ectoparasite sampling and identification 510
The present study identified 21 novel genetic lineages, plus 7 new ectoparasite species 511 records, from 138 bats of 17 species, sampled across 18 sites and 2 years. This suggests a 512 diverse ectoparasite fauna in this previously unsurveyed region of Mexico, but it is also likely 513 to be an underestimate of the true diversity, due to constraints around sampling effort. Bat 514 sampling was limited to May-August and conducted at relatively accessible locations with 515 water sources to facilitate bat capture. While our sampling sites were chosen to be 516 representative of habitat types across Baja, increasing the spatial and temporal scope of 517 sampling would be likely to increase the number of ectoparasite species discovered. 518 Assessment of ectoparasite fauna found in roosting sites against those found feeding directly 519 from their hosts and expanding seasonal coverage will be important for future work.