Climate change models predict southerly shift of the cat flea (Ctenocephalides felis) distribution in Australia

Background Bioclimatic variables play an integral part in the life-cycle of Ctenocephalides felis, the most common flea found on companion animals. It is essential that we understand the effects of climate on C. felis distribution as fleas are a major veterinary and public health concern. This study investigated the current distribution of C. felis in Australia and future projections based on climate modelling. Results Typing of C. felis was undertaken using the cytochrome c oxidase subunit 1 (cox1) mitochondrial DNA (mtDNA) region and current distribution of haplotypes was mapped by Maximum Entropy (Maxent) niche modelling. All C. felis haplotypes have been predicted to persist in environments along the eastern and southern coastlines of Australia and distinct ecological niches were observed for two C. felis haplogroups. Clade ‘Cairns’ haplogroup thrives under the northern coastal tropical conditions whilst Clade ‘Sydney’ haplogroup persists in temperate climates along the eastern and southern coasts. The model was then used to predict areas that are projected to have suitable climatic conditions for these haplogroups in 2050 and 2070 under the Intergovernmental Panel on Climate Change (IPCC) climate change scenarios. Under all IPCC Representative Concentration Pathways (RCP) climate change scenarios, the geographical range of all haplotypes was reduced by 5.59–42.21% in 2050 and 27.08–58.82% by 2070. The ranges of all clades were predicted to shift south along the eastern coastline. Conclusions As future temperatures exceed critical threshold temperatures for C. felis development in the northern tropical areas, Clade ‘Cairns’ haplogroup is predicted to shift south along the coastline and possibly outcompete the temperate haplogroup in these areas. If C. felis haplogroups possess distinct climatic niches it suggests a potential for these to be biologically distinct and have differing developmental rates and vector capabilities. Electronic supplementary material The online version of this article (10.1186/s13071-019-3399-6) contains supplementary material, which is available to authorized users.


Background
Ctenocephalides felis (Siphonaptera: Pulicidae), commonly known as the cat flea, is the most common flea found on companion animals in Australia [1,2]. Ctenocephalides felis has a cosmopolitan distribution and is highly tolerant to a wide range of environmental conditions [3]. As C. felis imposes health risks to humans and domestic animals as a biological vector, it is important to understand and predict both current and future suitable habitats [3]. Climatic variables play an integral part in the life-cycle of C. felis and consequently have an effect on their distribution, as they will only reside and successfully reproduce within limited climatic ranges [4]. As a result, it is essential that we understand the effects of climate on C. felis distribution, both now and in the future.
Recent studies have investigated the molecular diversity of C. felis in Australia at the cytochrome c oxidase subunit 1 (cox1) and cytochrome c oxidase subunit 2 (cox2) mitochondrial DNA (mtDNA) regions [1,2,[5][6][7]. Two genetically distinct subpopulations were identified, with one population located along the eastern and southern coasts and the other strictly located in the northern

Open Access
Parasites & Vectors *Correspondence: jan.slapeta@sydney.edu.au Sydney School of Veterinary Science, Faculty of Science, University of Sydney, Sydney, New South Wales 2006, Australia city of Cairns, Australia. It is currently unknown what biological factors are governing these distributions and whether the Cairns subpopulation is restricted to this area.
Climate change is anticipated to alter the ranges of zoonotic parasite species and consequently, associated disease emergence within human and domestic animal populations in the future [4,8]. The effect of climate change can have the potential to shift tropical populations into temperate areas as bioclimatic norms exceed critical threshold temperatures for parasite survival [9]. In a predictive model study, the distribution of C. felis in Spain was predicted to expand to newly suitable habitats as a result of climate change [10]. This could be the circumstance in Australia, where the distribution of C. felis subpopulations may shift to newly suitable habitats. Ecological niche modelling such as the Maximum Entropy (Maxent) model [8,11] in epidemiology is a useful tool as it can assess the relative importance of bioclimatic variables and use these factors to predict changes in the distribution of parasites and their pathogens over time [12].
The aim of this study is to model if C. felis distribution in Australia will be affected by the Intergovernmental Panel on Climate Change (IPCC) climate change scenarios. To address this aim, we evaluated the genetic diversity of the cat flea by polymerase chain reaction (PCR) amplification and sequencing of the cox1 mtDNA region. The predictive Maxent modelling was then used to determine the current and future distribution of C. felis haplotypes in Australia using IPCC climate change scenarios.

Flea specimens
Fleas were collected opportunistically from domestic cats and dogs presented to various veterinary clinics across the north-eastern region of Queensland (Additional file 1: Table S1). Fleas from individual animals were collected by veterinarians in 1.5 ml Eppendorf tubes, labelled with the postcode of the veterinary clinic, and stored in 80-100% ethanol at -20 °C for transport to The University of Sydney Veterinary Parasitology Diagnostic Laboratory. In addition, 65 C. felis samples from previous studies that have been published on GenBank and 33 unpublished C. felis genotyped samples were included (Additional file 1: Table S2) [1,2,5]. These additional samples were collected from veterinary clinics in all seven states and territories across Australia and were used to increase the sensitivity of the ecological niche model. There is a limitation in the extent of flea locality in this study as samples were collected from dogs or cats presented at veterinary clinics within their local area where these animals have acquired these fleas from an unknown locality.

Morphological identification of flea specimens
Fleas were identified morphologically using a stereo microscope (5-20× objective) with the aid of identification keys and descriptions [13,14].

Extraction of total DNA and mounting of flea exoskeletons
Total flea DNA was extracted using an ISOLATE II Genomic DNA kit (Bioline, Eveleigh, Australia) according to manufacturer's protocol with a few modifications as previously described [2]. Briefly, an incision was made on the anterior dorsal section of the flea abdomen using a sterile scalpel blade and digested overnight with 25 µl Proteinase K and 180 µl of lysis buffer at 56 °C. The remaining steps were as per manufacturers protocols. The total DNA was eluted into 50 µl of elution buffer and stored at -20 °C until analysis by PCR. Exoskeletons were placed in 10% potassium hydroxide for 12 h, rinsed with RO water and dehydrated in an ethanol series (70%, 80%, 95% and 100%) for 1 h each. The specimens were then slide-mounted in Euparal (Australian Entomological Supplies, Coorabell Australia).

Amplification of the cox1 gene by PCR and sequence analysis
The 601-nt fragment of the 5′-end of the cox1 gene sequence from the mtDNA region was amplified by PCR using a generic invertebrate amplification forward primer, LCO1490 [S0867] (5′-GGT CAA CAA ATC ATA AAG ATA TTG G-3′) and a previouslydesigned reverse primer, Cff-R [S0368] (5′-GAA GGG TCA AAG AAT GAT GT-3′) [2,15]. Each reaction was conducted at a final volume of 30 µl, containing 15 µl MyTaq ™ Red Mix (Bioline, Eveleigh, Australia) and approximately 2 µl of genomic DNA template. The PCR was conducted using C100 Thermal Cycler (BioRad, Gladesville, Australia) with cycling parameters set as following: denaturing at 95 °C for 1 min followed by 35 cycles at 95 °C for 15 s, 55 °C for 15 s and 72 °C for 10 s, followed by a final elongation for 5 min at 72 °C. All reactions were run with a negative control using sterile PCR-grade water [1,2]. Resulting amplicons were visualised on a 1% (w/v) agarose gel with GelRed ™ Nucleic Acid Gel Stain (Biotium, Fremont, CA, USA) in 1× TBE buffer to verify product size. Samples found to yield an unambiguous single band of expected size were sent for bidirectional sequencing (Macrogen Ltd, Seoul, Korea). Individual sequences were assembled and aligned against a reference sequence previously submitted to GenBank (KF684884) [2] using CLC Main Workbench 6.9.1 (CLC bio, Qiagen, Aarhus, Denmark).

Modelling of current C. felis distribution in Australia using Maxent
A Maxent model [8,11] was used to model the current distribution of C. felis in Australia. A total of 179 samples were used, which included samples obtained in this study, unpublished data and published data (Additional file 2: Table S4). Latitude and longitude coordinates for each sample were entered into Biodiversity and Climate Change Virtual Laboratory (BCCVL) (BCCVL, Griffith University) as biological data [16]. "Australia, Current Climate (1976-2005), 30 arcsec (~1 km)" climate and environmental data was selected from the datasets available in the BCCVL collections. All variables were initially included in the model, while only the variables with the highest probability of C. felis presence in response to ecogeographical variables were included in the final model as they reported the greatest amount of the observed variation. The chosen response variables were required to be within the environmental thresholds for C. felis survival. Namely, temperatures between 3-35 °C with 70-95% humidity and high precipitation levels (> 500 mm yearly) [17,18]. Receiver operating characteristics (ROC) curves were assessed, where the area under the curve (AUC) was used to evaluate the accuracy of the resulting model. A model with an AUC score of 0.5 or below indicated that the model performed no better than random. Models with an AUC score of 1 indicated a perfect model. Additionally, the variable importance function was used to calculate the correlation score between the standard prediction and the new prediction, giving an estimation (%) of the importance of that variable in the model. The higher the value of the variable, the more influential it was in the model. The correlation matrix was used to complement the variable importance function where it showed the probability of predictor variables to be correlated to ensure that the included variables are not highly correlated and produce a bias output in the model. It is important to note that there are limitations involved in strictly using AUC models, and they should be used in conjunction with other model evaluation methods [19]. Therefore, we used the variable importance function and correlation matrix with the AUC scores to determine the model of best-fit for the current distribution of C. felis.

Species distribution modelling confirms a geographical niche for Clade 'Cairns' and the nation-wide distribution of Clade 'Sydney'
An additional 98 samples were used in conjunction with the 81 samples obtained from this study to model the species distribution of C. felis in Australia (n = 179; Fig. 1). The best Maxent model included a fine scale with the following four environmental variables: maximum temperature of warmest month, mean temperature of coldest quarter, precipitation of wettest quarter and precipitation of warmest quarter (Additional file 3: Figure S1). The best model was determined by interpreting AUC values, variable importance functions and the correlation matrix ( Fig. 2). The AUC for Clade 'Sydney' , Clade 'Darwin' , Clade 'Cairns' and all haplotypes were 0.96, 0.99, 1.0 and 0.98, respectively. The variable importance function identified that the maximum temperature of warmest month, mean temperature of coldest quarter, precipitation of wettest quarter and precipitation of warmest quarter variables had an importance of 10%, 40%, 25% and 25%, respectively, in the model. The correlation matrix showed that precipitation of wettest quarter and precipitation of warmest quarter are highly correlated (0.5) whereas all other variables are uncorrelated with each other.
Areas that had a suitability > 50% for Clade 'Sydney' , Clade 'Darwin' and Clade 'Cairns' had a maximum temperature of 35 °C, 40 °C and 35 °C, respectively, in the warmest month, mean temperature of 10 °C, 15 °C and 25 °C, respectively, in the coldest quarter, precipitation of 250 mm, 1700 mm and 500 mm, respectively, during the wettest quarter and lastly, a precipitation of 1400 mm, 1400 mm and 1500 mm, respectively, during the warmest quarter (Fig. 2). These conditions fall within the required thresholds for survival of C. felis. Clade 'Sydney' was observed to be found along the entirety of east coast of Australia, extending to the coast of South Australia, Tasmania and Perth (Fig. 2). Areas with suitability for the Clade 'Darwin' were primarily found along the northern coastline of Australia, the north to middle parts of the east coast and the western coast of Tasmania (Fig. 2). In contrast, suitable habitats for Clade 'Cairns' were primarily found in the northern and north-eastern tropical areas of Australia, in particular Townsville, Cairns and the surrounding region (Fig. 2).

Modelling of potential future distribution predicts the range of C. felis in Australia to shift south along the eastern coastal regions
Areas with suitable climatic conditions for C. felis in Australia have been predicted to decrease by 2050 and even further by 2070 under all four RCP scenarios (2.6, 4.5, 6.0 and 8.5; Table 1). Currently, 959,040.82 km 2 is suitable for C. felis in Australia. However, in 2050 the suitable area reduced in all four scenarios, ranging between 548,751.36-680,130.58 km 2 , leaving 54.45-68.61% of the area in common with the current model (Table 1, Fig. 3). By the 2070s, the area of suitability declined even further in all four models with 386,233.19-717,111.65 km 2 suitable for C. felis, leaving 27.82-70.26% of the area in common with the current model (Table 1, Fig. 3). The contracted areas are found in the current northern parts of the C. felis range where the geometric centre of the species range shifted 397.14 km south along the east coast by the 2050s (Table 1; Fig. 3). This southward shift is observed to be a further 127.47 km by the 2070s (Table 1, Fig. 3 Fig. 4). This resulted in a subsequent shift rate of 122.11 km south per decade, with loss of the current northern range (  Fig. 6).

Discussion
This study has doubled the number of fleas genotyped in Australia and provided an updated perspective on the diversity of C. felis, particularly in the north-eastern region. Only three cox1 haplotypes have previously been reported from 47 C. felis specimens collected from domestic cats and dogs across Australia [1,2,5,7]. This study confirms the presence of the haplotypes (h1, h4 and h5) previously genotyped, and uncovers the existence of three new cox1 haplotypes (h2, h3 and h6) in north-eastern Australia. These results provide strong evidence supporting intra-specific diversification of C. felis between  the north-eastern region and middle to southern regions of Australia. As haplotypes h4, h5 and h6 have only been detected in this tropical region and differ genetically by one SNP, they have been categorised as Clade 'Cairns' , whilst genetically similar haplotypes h1 and h2 have been found elsewhere have been grouped as Clade 'Sydney' . The findings from this study suggest a range expansion of tropical haplogroups as only Clade 'Sydney' haplogroup was previously observed in the surrounding area of Cairns in 2014 [1,2]. However, the theory of emerging haplotypes is questionable as the one dog that was sampled in the outer area of Cairns may not have presented the emerging haplotypes [1,2]. It was unexpected to observe haplotype h3 in the Cairns region as previously it has only been found in Galiwinku, near Darwin (our unpublished data). Clade 'Darwin' may have been present throughout the northern tropical regions before this study and not sampled as only 25 samples have been previously genotyped in these regions.
An initial bottleneck restricting gene flow of the flea populations into Australia could have been the cause of the low haplotype diversity currently present [24]. Over time, the result of extensive transportation and the movement of cats and dogs has facilitated interaction between flea populations [25]. Additionally, the abundance of feral cats and dogs could expedite greater genetic transfer between flea populations [25]. Despite an increase in diversity suggested by this study, it still reveals relatively low mitochondrial variability and highly restricted haplogroup populations of C. felis in Australia. The finding of six mtDNA haplotypes is considered to be unusually low compared to other flea species, such as the ten haplotypes of Pulex simulans and 24 haplotypes of Oropsylla hirsuta, the prairie dog flea [26,27]. It is expected that for at least every ten samples of mtDNA sequences greater than 500 bp (e.g. cox1), multiple haplotypes are expected for one species of animal [28]. Nevertheless, low mitochondrial variability has also been observed in other arthropod species [29]. High homology was similarly found at the cox1 and cox2 mtDNA regions in a study of Anopheles stephensi stephensi, an important malarial vector in India [29]. Oshaghi et al. [29] suggested the low genetic variation is a result of unrestricted gene flow due to the lack of geographical barriers. This, along with extensive movement of hosts, could be the potential contributing factor for the low genetic variation found within Australia [30].
The extensive movement of host allows indirect dispersal of fleas over a large range; however, the results from this study and two previous studies [1,2] suggest three genetically distinct haplogroups of C. felis in Australia. These genetically distinct groups indicate a lack of gene flow between some populations [31]. The species diversification is likely to be influenced by climatic variation as climatic differences can impede gene flow despite sufficient opportunity for dispersal from frequent movement of hosts [4].
The model for the present distribution of C. felis in Australia shows that the most suitable habitat is located along the coastline as it presents the ideal environment for flea survival. The environmental conditions of the eastern coast of Australia consists of moderate to warm temperatures and high levels of precipitation that are within the critical thresholds for C. felis survival (Australian Bureau of Meteorology, Australian Government; http://www. bom.gov.au/iwk/clima te_zones /map_1.shtml ) [17,18]. The Clade 'Darwin' model revealed that the range of haplotype h3 is along the northern coastal environments. This model for Clade 'Darwin' requires additional samples from these remote and scarcely populated areas to be included to further refine the model. Ctenocephalides felis does not have the ability survive and develop in temperatures greater than 35 °C as is seen in inland regions [32].
The predictive models indicate that there is a tropical ecological niche for Clade 'Cairns' and Clade 'Darwin' , whereas Clade 'Sydney' has the ability to inhabit most Australian climates. Species diversification as a result of climate as a natural barrier has been observed in other species. The climatic variation between the northern and southern part of Fukushima and Ibaraki prefectures in Japan, with mean temperatures of 18.3 °C and 17 °C respectively, is thought to be the reason for the presence of two genetically distinct groups of the sorghum plant bug, Stenotus rubrovittatus [33]. Australia has a very variable climate across the continent where the northern region presents a tropical climate that transitions into a temperate climate in the southern states (Australian Bureau of Meteorology, Australian Government; http:// www.bom.gov.au/). Fleas dispersed from host movement may not survive, reproduce or parasitise hosts outside their origin due to movement to a new unfavourable habitat [34].
As climate is a fundamental driver in shaping species diversification, climate change can lead to a range shift or reduction in habitat suitability [35]. Under all RCP climatic scenarios, a southward shift in the distribution of C. felis was observed to occur along the coastal areas of Australia by the 2050s and continuing further by the 2070s (range shift = 124.11 km per decade). Although the proximate expansion of C. felis into new suitable regions continues by extensive movement of their hosts, our findings predict an ultimate decline of 61.43% in suitable climatic habitat by the 2070s under the RCP8.5 model. The decrease in suitable habitat further north is most likely due to a predicted increase in the maximum temperature of the warmest month. An increase in temperature can increase the rate of development for arthropods; however, the positive association between temperature and developmental rate may be offset by the species total bioclimatic requirement for survival [36]. Initially, the increased temperatures may promote the spread and occurrence of C. felis along the north-eastern coastline, but the expected temperature increase of 1.1-6.4 °C from 1980 to 2099 will exceed critical threshold temperatures (35 °C) for C. felis survival in the current northern areas of habitat [37]. As a result, the climatic requirements for the growth and development of C. felis will become limited in these areas and consequently restrict the distribution of the parasite to areas further south with sufficient bioclimatic resources [37]. Climate change in the future is also expected to increase precipitation and humidity levels in some areas while other areas will experience severe drought conditions [38]. It is known that flea infestations are mostly absent in the Australian inland communities as drought conditions are too harsh for flea survival [39]. As C. felis development and survival is highly dependent on moist environments, distribution is most likely to be found in those regions [17]. As moist regions are becoming further restricted to coastlines in the future, the model suggests that there will be an increase in flea populations in these areas.