Spatial Distribution of Forensically Significant Blow Flies in Subfamily Luciliinae (Diptera: Calliphoridae), Chiang Mai Province, Northern Thailand: Observations and Modeling Using GIS

Blow flies of the subfamily Luciliinae (Diptera: Calliphoridae) are one of the main forensically important subfamilies globally. In addition to being used to estimate the minimum post-mortem interval (PMImin), assuming colonization occurred after death, blow fly specimens found infesting a human corpse are used to determine if the corpse was relocated or if the individual ingested narcotics prior to death. The presence of these blow flies in a given area is strongly influenced by abiotic and biotic factors, such as temperature, elevation, and habitat. Having this information, along with geographical distributions and the characteristics of preferred habitats, is necessary to better understand the biology of this group. This study aimed to characterize the spatial distribution of Luciliinae throughout 18 sampling sites within six ecozones (disturbed mixed deciduous forest, mixed deciduous forest, mixed orchard, paddy field, lowland village, and city/town) in central Chiang Mai Province, northern Thailand over one year (May 2009–May 2010). The purpose of the study was to elucidate the relationship of blow fly species composition with environmental abiotic factors (e.g., temperature, relative humidity, light intensity), and to predict the distribution of the common species within this subfamily using GIS. Adult collections were performed biweekly, baited with one-day-old beef offal. A total of 2331 Luciliinae flies trapped, comprising eight species, of which the four predominant species were Hemipyrellia ligurriens (Wiedemann) (n = 1428; 61.3%), Lucilia porphyrina (Walker) (n = 381; 16.3%), Hemipyrellia pulchra (Wiedemann) (n = 293; 12.6%), and Lucilia papuensis Macquart (n = 129; 5.5%). Population density across species varied seasonally, peaking in August 2009 coinciding with the rainy season. Predicting population composition was based on a model developed with ArcGIS 9.2, which utilized environmental variables (temperature, relative humidity, and light intensity) in conjunction with abundance data. Models indicated H. ligurriens had the most widespread geographic distribution, while H. pulchra was predicted to occur largely in mixed orchards and lowland villages. Lucilia porphyrina and L. papuensis were less widespread, restricted mainly to mixed deciduous forest. This model, along with knowledge of forensic information, may be useful under certain investigations where the corpse may have been relocated.


Adult Flies Collection
Adult flies were trapped every two weeks from May 2009 to May 2010, using an in-house prototype portable funnel trap kit [33]. The trapping method was described previously [33]. Briefly, the trap consisted of a polyvinyl chloride (PVC) frame box (30 × 30 × 50 cm), a fly entrance module, and a black fly net (30 × 30 × 80 cm). Two replicate sets of traps were placed in a single row (50 m apart) on the ground [33], each baited with 250 g of one-day-old beef offal [37], which was kept in a translucent plastic container. The bait was placed underneath the fly entrance module after it had been tied to a fly net with elastic bands and installed with the PVC frame box. Flies lured by the smell from the offal were collected by their passive movement up toward the light through the fly entrance module to the fly collection net after landing on the bait. Traps were placed in the shade to prevent thermal stress for trapped flies [38]. The collection was done for a 1-h period between 09:30 a.m. and 12:00 p.m. Physical data of each study site were collected including light intensity (lux) (LUX/FC light meter TM-204 Tenmars, Tenmars Electronics Co., Ltd., Taipei, Taiwan), temperature (°C), relative humidity, RH (%) (Digital Hygro-Thermo Meter (DHT-1), Daeyoon Scale Industrial Co., Ltd., Seoul, South Korea) and the co-ordinates (Garmin™ eTrex Handheld GPS, Garmin China Co., Ltd., Chaoyang, China).
All specimens were carried to the laboratory of the Department of Parasitology, Faculty of Medicine, Chiang Mai University. They were sacrificed by placement in a freezer set at 0 °C for 2 h. Flies were then identified under a dissecting microscope (model SZ2-ILST, Olympus Corporation, Tokyo, Japan) using taxonomic keys, sexed, and counted [20,21].

Statistical Analysis
Three seasons were considered: a rainy season occurring from June through October, winter from November through February, and summer from March through May (File S2). Species sampled with fewer than 100 specimens captured year-round were excluded from the analysis. We analyzed

Adult Flies Collection
Adult flies were trapped every two weeks from May 2009 to May 2010, using an in-house prototype portable funnel trap kit [33]. The trapping method was described previously [33]. Briefly, the trap consisted of a polyvinyl chloride (PVC) frame box (30 × 30 × 50 cm), a fly entrance module, and a black fly net (30 × 30 × 80 cm). Two replicate sets of traps were placed in a single row (50 m apart) on the ground [33], each baited with 250 g of one-day-old beef offal [37], which was kept in a translucent plastic container. The bait was placed underneath the fly entrance module after it had been tied to a fly net with elastic bands and installed with the PVC frame box. Flies lured by the smell from the offal were collected by their passive movement up toward the light through the fly entrance module to the fly collection net after landing on the bait. Traps were placed in the shade to prevent thermal stress for trapped flies [38]. The collection was done for a 1-h period between 09:30 a.m. and 12:00 p.m. Physical data of each study site were collected including light intensity (lux) (LUX/FC light meter TM-204 Tenmars, Tenmars Electronics Co., Ltd., Taipei, Taiwan), temperature ( • C), relative humidity, RH (%) (Digital Hygro-Thermo Meter (DHT-1), Daeyoon Scale Industrial Co., Ltd., Seoul, South Korea) and the co-ordinates (Garmin™ eTrex Handheld GPS, Garmin China Co., Ltd., Chaoyang, China).
All specimens were carried to the laboratory of the Department of Parasitology, Faculty of Medicine, Chiang Mai University. They were sacrificed by placement in a freezer set at 0 • C for 2 h. Flies were then identified under a dissecting microscope (model SZ2-ILST, Olympus Corporation, Tokyo, Japan) using taxonomic keys, sexed, and counted [20,21].

Statistical Analysis
Three seasons were considered: a rainy season occurring from June through October, winter from November through February, and summer from March through May (File S2). Species sampled with fewer than 100 specimens captured year-round were excluded from the analysis. We analyzed the relationship of environmental variables (temperature, relative humidity, and light intensity) with trap Insects 2018, 9, 181 4 of 15 catch using negative binomial regression analysis. We used one-way analysis of variance (ANOVA) followed by the post hoc Bonferroni test (homogeneity of variance: p > 0.05) or Dunnett T3 test (homogeneity of variance: p < 0.05) to determine the relationship between land use types and mean total number of flies using IBM SPSS Statistics for Windows, version 22.0 (IBM Corp., Armonk, NY, USA) (α = 0.05).
To predict the spatial distribution of Luciliinae flies by season and aggregated for the full year, we conducted additional geospatial analysis using a kriging/co-kriging approach in ArcGIS 9.2. Kriging is a kind of linear least squares estimation method [39] that is associated with spatial optimal linear prediction in which the unknown random-process mean is estimated with the best linear unbiased estimator [40]. To create a continuous surface of the phenomenon, predictions are made for each location in the study area based on the semivariogram and the spatial arrangement of measured values that are nearby. Co-kriging is an extension of kriging for prediction of one variable using other variables. The co-variables must have a relationship and this relationship must be defined [39]. The use of co-kriging requires the spatial covariance model of each variable and the cross-covariance model of the variables is defined [39].
In this study, the parameters that were used in co-kriging included climatic factors that had a statistically significant relationship with fly numbers and the total trap catch during the study period. The data were log-transformed to normalize the distribution and minimize standard error of the geographical analysis [40]. Moreover, the data of land use types were defined as dummy variables and were incorporated into the ordinary kriging/co-kriging analysis. The mathematical models for evaluating the semivariogram/covariance function were accepted when the model provided the lowest root-mean-square prediction error [33].

Results
A total of 2331 Luciliinae flies were sampled, comprising eight species, of which the predominant species was H. ligurriens, representing 61.3% of the individuals (n = 1428). Other associated species were L. porphyrina (n = 381; 16.3%), H. pulchra (n = 293; 12.6%), and L. papuensis (n = 129; 5.5%). Small numbers of L. cuprina (n = 67; 2.9%), Hypopygiopsis tumrasvini Kurahashi (n = 23; 1.0%), Lucilia sinensis Aubertin (n = 6; 0.3%) and Hypopygiopsis infumata (Bigot) (n = 4; 0.2%) ( Table 1). Therefore, this paper focused mainly on the four most sampled species, namely H. ligurriens, L. porphyrina, H. pulchra, and L. papuensis. The sampled fly numbers varied seasonally, peaking in August 2009, coinciding with the rainy season, while a minor peak was in summer (May 2009). The sample numbers decreased gradually during the late rainy season (September 2009) and remained low throughout the winter and summer ( Figure 2).     Hemipyrellia ligurriens had widespread geographic distribution (Table 1). Although this species was sampled throughout the year, a bimodal peak was observed-the major peak occurring in late summer (May 2009), while a minor peak occurred in the rainy season (June-October 2009). As expected, sample numbers gradually decreased at the onset of winter (November 2009) and remained minimal throughout winter and summer (January-May 2010) ( Figure 2). H. ligurriens trap catch was positively correlated with temperature and relative humidity (p < 0.05), with increased temperature (B = 0.103, p = 0.0001) and relative humidity (B = 0.035, p = 0.0001) associated with increases in trap catch ( Table 2). No correlation with light intensity was found (p = 0.954). The greatest numbers were collected at 25-30 • C and 60-70% RH (Figure 3). Analyses using co-kriging showed the spatial distribution of fly density across a variety of land uses in Mueang Chiang Mai district ( Figure 4).     Lucilia porphyrina was prevalent in mixed deciduous forest (Table 1)  Lucilia porphyrina was prevalent in mixed deciduous forest (Table 1) Hemipyrellia pulchra was collected primarily in mixed orchard and lowland village; this species was not found in paddy field or city/town (Table 1) Table 2). The greatest number of this species was captured at 25-30 °C and 60-70% RH (Figure 3). Co-kriged prediction maps were produced, showing abundance in mixed orchard and lowland village. In summer, the spatial pattern was similar to the aggregated one-year survey. This species was predicted to occur largely in mixed orchard and disturbed mixed deciduous forest in the rainy season, while the fewest numbers were predicted in winter, with little apparent spatial variation ( Figure 6). Hemipyrellia pulchra was collected primarily in mixed orchard and lowland village; this species was not found in paddy field or city/town (Table 1) Table 2). The greatest number of this species was captured at 25-30 • C and 60-70% RH (Figure 3). Co-kriged prediction maps were produced, showing abundance in mixed orchard and lowland village. In summer, the spatial pattern was similar to the aggregated one-year survey. This species was predicted to occur largely in mixed orchard and disturbed mixed deciduous forest in the rainy season, while the fewest numbers were predicted in winter, with little apparent spatial variation ( Figure 6). Lucilia papuensis was significantly prevalent (ANOVA, p < 0.05) in mixed deciduous forest and lowland village, respectively (Table 1) Table 2). The greatest numbers of flies were collected at 25-30 °C and 70-80% RH (Figure 3). Co-kriged prediction maps were produced. A one-year survey and seasonal prediction maps exhibited a similar pattern, as this species was predicted to occur largely in mixed deciduous forest along the mountainous areas and in lowland villages (Figure 7). Lucilia papuensis was significantly prevalent (ANOVA, p < 0.05) in mixed deciduous forest and lowland village, respectively (Table 1). A bimodal peak population was evident, with a major peak in the rainy season (August-September 2009) and a minor peak at the onset of the rainy season  Table 2). The greatest numbers of flies were collected at 25-30 • C and 70-80% RH (Figure 3). Co-kriged prediction maps were produced. A one-year survey and seasonal prediction maps exhibited a similar pattern, as this species was predicted to occur largely in mixed deciduous forest along the mountainous areas and in lowland villages (Figure 7).

Discussion
An adequate knowledge of the ecology and geographical abundance of blow flies is important for the forensic sciences [2,3,41]. However, in Thailand, such information on geographic distribution of forensically important flies was limited prior to this study, which is the first to report the distribution of four forensically important Luciliinae flies; H. ligurriens, H. pulchra, L. porphyrina, and L. papuensis, over a full year. Maps of predicted fly distribution were generated using the Geostatistical Analyst Extension of ArcGIS 9.2.
Abundance relative to ecotype differed among the four species. The study sites, which represented different broad ecological locations, were chosen based on a systematic random sampling method. Some of the 18 study sites under the six land use categories, covering the city to forested areas, may not be suitable for some fly species [17,20,21,42]. We found H. ligurriens had widespread geographic distribution, while H. pulchra was collected mainly in mixed orchard and lowland village land uses. This result may reflect that H. ligurriens can live and develop in many substances such as animal dung, carcasses, human feces, and decomposed organic matter [17,19], while the adult of H. pulchra prefer flowers and fruits in the orchards [19]. A study in Phitsanulok Province, northern Thailand, also reported a wide variety of habitat types, where H. ligurriens can be found, ranging from residential areas to the mountainous zones, but H. pulchra was collected only in

Discussion
An adequate knowledge of the ecology and geographical abundance of blow flies is important for the forensic sciences [2,3,41]. However, in Thailand, such information on geographic distribution of forensically important flies was limited prior to this study, which is the first to report the distribution of four forensically important Luciliinae flies; H. ligurriens, H. pulchra, L. porphyrina, and L. papuensis, over a full year. Maps of predicted fly distribution were generated using the Geostatistical Analyst Extension of ArcGIS 9.2.
Abundance relative to ecotype differed among the four species. The study sites, which represented different broad ecological locations, were chosen based on a systematic random sampling method. Some of the 18 study sites under the six land use categories, covering the city to forested areas, may not be suitable for some fly species [17,20,21,42]. We found H. ligurriens had widespread geographic distribution, while H. pulchra was collected mainly in mixed orchard and lowland village land uses. This result may reflect that H. ligurriens can live and develop in many substances such as animal dung, carcasses, human feces, and decomposed organic matter [17,19], while the adult of H. pulchra prefer flowers and fruits in the orchards [19]. A study in Phitsanulok Province, northern Thailand, also reported a wide variety of habitat types, where H. ligurriens can be found, ranging from residential areas to the mountainous zones, but H. pulchra was collected only in the agricultural areas and in the forests [43]. A previous study reported various altitude distribution ranges of H. ligurriens, indicating greater adaptation compared with H. pulchra [32].
Lucilia porphyrina and L. papuensis were less widespread, and were restricted mainly to the mixed deciduous forest (at 950 m a.s.l.). Likewise, others have reported that L. porphyrina and L. papuensis are likely to have their distribution restricted to forested highland areas in Thailand [21,32,43]. A study in Japan indicated that the females of L. porphyrina preferred to lay their eggs in the mountains and/or highlands [44]. Additionally, their larvae have been reported to feed on animals' carcasses [19,21], which could be found in the forests. However, L. porphyrina could be found in human residences in India and in the mainland of Japan (for overwintering) [19,44]. A small number of L. papuensis was collected in the present study, but it was found in all six of the land uses, with the greatest number in mixed deciduous forest. Many studies also supported its preference for the forest areas [17,19,42,43]. This may be caused by the adults gathering around decaying animals (e.g., earthworms, land snails, and snakes) and other vertebrate carcasses [19,21] that could be found in the forest areas.
The peak trap catch of these species was found in the rainy season (August 2009), and the nadir occurred during the summer (March-May 2010). This finding is in contrast with many studies where more blow flies were collected during the warmer months [22,45,46]. The conditions in the rainy season (mean temperature 29.1 ± 0.3 • C, mean RH 69.7 ± 1.0%) may be suitable for the development of Luciliinae flies. A previous study with H. ligurriens suggested the optimal condition for the development of this species is a temperature range of 16 to 28 • C, not 31 • C or above [47]; therefore, the temperature of the summer months in this study (mean temperature 30.6 ± 0.5 • C) (see Figure 2) may exceed the optimal range. However, the present study was carried out during a single year. To obtain reliable seasonal patterns, the season needs to be replicated (i.e., the study should be replicated during more than one year). In this study, fly trapping was done before noon (9:30 a.m. and 12:00 p.m.). Blow fly activity was reported to be higher around noontime and in the afternoon than in the early morning [22,48,49], so our collection time may have affected the resulting collections.
The yearly sampling across 18 study sites allows us to spatially predict the population distribution of Luciliinae flies. Our study is not the first to model the prediction of forensically important blow flies. Previous investigations have shown the spatial models of C. megacephala, C. rufifacies, C. pinguis, C. chani, C. villeneuvi and Cey. nigripes [33][34][35]. According to this study, H. ligurriens was predicted to occur across a variety of land uses in Mueang Chiang Mai district especially nearby residential areas (e.g., city/town and lowland village), which is similar to C. megacephala and C. rufifacies [33,34]. On the other hand, L. porphyrina and L. papuensis showed preference in the mixed deciduous forest, which is consistent with C. pinguis, C. chani, C. villeneuvi and Cey. nigripes [35]. The greatest number of H. pulchra was predicted mainly in mixed orchard and lowland village, with low altitudes ranging between 300 and 400 m a.s.l., which is in contrast with a previous study that found this species only at 901-1050 m a.s.l. [32]. However, their studies were carried out in the forested area of altitudes ranging from 400-1050 m a.s.l. and only one specimen of H. pulchra was collected by a sweep net.
According to our results, the prediction maps indicate H. ligurriens occurred in a variety of habitats and is closely associated with human activity. Previous studies reported H. ligurriens infested a human corpse in a forested area of northern Thailand [4]. This species has also been associated with human remains in an outdoor environment of Australia [11]. The widespread distribution of H. ligurriens (ranging from residential areas, forests, and highland) was also observed in Bangladesh and Malaysia [7,42]. In contrast, L. porphyrina was predicted mainly in the mixed deciduous forest at 905 m a.s.l. This species may be used as insect evidence when corpses are found in the forest of mountainous areas. Recently, L. porphyrina was found on a human corpse discovered in an open mountainous area (1200 m a.s.l.) during the winter time of northern Thailand [16]. Additionally, this species infested rabbit carcasses in the highlands of Malaysia (1517.3 m a.s.l.), suggesting that it is a forensically important species in the highlands.

Conclusions
We provide predictive maps of the spatial distribution for forensically important blow flies within the subfamily Luciliinae as related to land use types and climatic factors. Hemipyrellia ligurriens and H. pulchra have spread through a wide range of land use types, implying that their application for explaining the location of corpse movement is limited, especially when the immature specimens of these species were found. On the other hand, L. porphyrina and L. papuensis evidently show a preference to the mixed deciduous forest. This database may be incorporated as a useful tool in forensic entomology. To be used for PMI min estimation, biological information such as developmental rate or insect succession of these species are mandatory.