國立臺灣大學生物資源暨農學院森林環境暨資源學系 碩士論文
School of Forestry and Resource Conservation College of Bioresources and Agriculture
National Taiwan University Master Thesis
使用遙測資料探討都會地景蜘蛛多樣性之空間分布 Spatial Patterns of Spider Diversity in Urban
Landscape by Using Remote Sensing Derived Indices
尤光平 Kuang-Ping Yu
指導教授：丁宗蘇 博士 及 卓逸民 博士 Advisors: Tzung-Su Ding, Ph.D. & I-Min Tso, Ph.D.
中華民國 108 年 6 月
June, 2019
i
國立臺灣大學碩士學位論文 口試委員會審定書
使用遙測資料探討都會地景蜘蛛多樣性之空間分布
Spatial Patterns of Spider Diversity in Urban Landscape by Using Remote Sensing Derived Indices
本論文係尤光平君（R06625007）在國立臺灣大學森林環
境暨資源研究所完成之碩士學位論文，於民國 108 年 6 月 5 日
承下列考試委員審查通過及口試及格，特此證明
ii
誌謝
首先，我要感謝丁宗蘇老師和卓逸民老師耐心的指導，從一開始研究題目的 選擇、後續的統計分析、到論文的撰寫修改與口試，大至研究、小至生活，兩位 老師都盡他們所能的在各個層面給予我非常多的協助，也在我兩年的碩士生涯中 給了我相當大的自由度。另外，也非常謝謝李培芬和張琪如兩位老師願意擔任我 的口試委員，除了指出我研究中的盲點之外，也給予了我很多實用的建議。
接著，謝謝研究室伙伴林穆明、呂立中、李先祐、蔡芷怡、林佳祈以及沈芳 伃在研究上提供的協助、建議，也謝謝張瀚柏、張晉嘉在文獻上的支援，希望大 家接下來的研究生涯都非常順利。另一方面，非常感謝東海大學黃博森學長每次 都不厭其煩的聆聽我提出關於原始資料的問題，也要謝謝菲菲學姐時常協助我尋
找卓老師並和老師聯絡，此外，謝謝葉浩然、吳婉瑜以及Sean Kelly 在研究以及
英語相關的問題給予我協助，沒有大家，這份論文不會如此順利的完成。
最後，希望好好謝謝我的家人，儘管我走上這條與普世價值相左的道路，並 一頭栽進了蜘蛛的世界，你們都還是全力支持我，並提供我各式各樣的支援，甚 至還願意在平日的晚上陪我去荒山野嶺採集蜘蛛，真的非常感謝你們。
iii
中文摘要
遙測技術可以監測大空間尺度下的環境變化，使探討都市化對生物多樣性的 衝擊研究更容易進行，其中，又以能夠反映棲地枝葉量多寡與初級生產力的植被 指數最常被使用於描述棲地的面積以及棲地的品質；然而，有關棲地初級生產力 假說和遙測資料衍生的環境因子如何影響蜘蛛這類肉食性小型節肢動物在大空間 尺度下分布的研究卻較為缺乏且結果分歧。因此，本研究使用遙測資料所衍生的 棲地動態指數與建蔽/裸地動態指數，探討與蜘蛛物種多樣性、科組成以及棲地偏
好間的關係，並找出遙測因子能解釋最多變異的取樣尺度。本研究透過1145 組
掉落式陷阱蒐集蜘蛛樣本，並將蜘蛛物種多樣性依照1、6.25、25、100 公頃的網
格進行資料整合；另一方面，遙測資料則由大地8 號衛星的網路公開資料庫取
得，並使用廣義線性混合模式以及冗餘分析測試遙測因子和蜘蛛物種多樣性資料 之間的關係。本研究發現，蜘蛛物種多樣性和棲地初級生產力呈顯著的正相關，
且由於上行效應，使得蜘蛛儘管是食物鏈中上層的掠食者而非直接取食植物的初 級消費者，棲地生產力假說仍能適用於牠們；且遙測因子，特別是年累積生產
力，能夠顯著地解釋蜘蛛的物種多樣性以及科組成，而解釋變異的比例在6.25 到
25 公頃的空間尺度之間達到最大。此外，狼蛛、皿蛛、姬蛛以及卵蛛這四個優勢 科基於其體型、覓食策略和型態特徵的不同，棲地偏好能良好地被遙測因子所解 釋。本研究提供了使用大尺度監測工具研究小型肉食節肢動物的方法，在保育策 略方面，於都會綠地中營造初級生產力高且穩定的環境有助於蜘蛛多樣性的保 育，除大面積綠地之外，高度都市化地區的行道樹或小型花壇亦是都會蜘蛛的重 要棲地，減少殺蟲劑和除草劑的用量對於都會蜘蛛的物種多樣性保育亦有助益。
關鍵字：棲地生產力、增強植被指數、棲地動態指數、增強建蔽/裸地指數、臺灣
iv
ABSTRACT
Through remote sensing data, monitoring the impact of urbanization on species
diversity at a regional scale has become more and more convenient. Vegetation indices
have been used to represent the size of habitats and habitat quality by predicting the
amount of foliage and habitat productivity. However, studies on how habitat primary
productivity hypothesis and remote sensing derived environmental factors (RS factors)
affect spatial distribution of spiders at large spatial scales are lacking, and the relationship
between spider diversity and RS factors remains unclear. Hence, the present study aimed
to examine the relationship of spider species richness, family composition, and habitat
preference with RS factors (Dynamic Habitat Indices and Dynamic Building/Bareness
Indices) and determine the best spatial scale of sampling unit which RS factors could
explain the largest variance in spider species richness and family composition. Spider
species distribution data were obtained by pitfall traps in 1,145 sampling sites in an
urbanization landscape in central Taiwan. Remote sensing data were obtained from
Landsat 8 images. The relationships between RS factors and spider assemblage diversity
were examined by generalized linear mixed models and redundancy analysis at four
spatial scales: 1, 6.25, 25, and 100 ha grids. Results reveal that although spiders are
predators which occupy higher trophic level, spider diversity follows habitat productivity
v
hypothesis based on bottom-up effect, and thus could be modeled by RS factors
significantly, especially cumulative Dynamic Habitat Index. The best spatial scale for
studying spider diversity by RS factors was between 6.25 to 25 ha. With differences in
sizes, foraging strategy, and morphological traits, habitat preferences of dominant spider
families, Lycosidae, Linyphiidae, Theridiidae, and Oonopidae, could also be well
explained by RS factors. Overall, the present study offers methods of modeling the spatial
distribution of small carnivorous invertebrate species richness in the urbanization
landscape by remote sensing data at a broad scale. For spider biodiversity conservation,
maintaining high and stable habitat productivity in green areas, such as parks or school
campuses, can help maintain high spider species richness. Also, Reducing the use of
pesticides and herbicide at street trees and small vegetated patches in highly urbanized
areas may also help on conserving spider diversity.
Keywords: Habitat productivity, EVI, Dynamic Habitat Indices, EBBI, Taiwan
vi
CONTENTS
口試委員審定書 ... i
誌謝 ... ii
中文摘要 ... iii
ABSTRACT ... iv
CONTENTS ... vi
LIST OF FIGURES ... vii
LIST OF TABLES ... x
Introduction ... 1
Methods ... 12
Results ... 24
Discussion ... 35
References ... 42
Appendix ... 54
vii
LIST OF FIGURES
Figure 1 The study area located in central Taiwan; each dot represent one sampling site and colors of the dots represent different land use types in Global Research &
Education on Environment and Society project. ... 13
Figure 2 Dynamic Habitat Indices (DHIs) of the study area derived by Enhance Vegetation Index. (a) cumulative DHI, (b) minimum DHI, (c) variation DHI,
(d) combined DHIs were transformed into RGB color on the map (Red:
Variation DHI; Green: Cumulative DHI; Blue: Minimum DHI)... 19
Figure 3 Dynamic Building/Bareness Indices (DBIs) of the study area derived by Enhance Build-up/Bareness Index. (a) cumulative DBI, (b) maximum DBI,
(c) variation DBI, (d) combined DBIs were transformed into RGB color on
the map (Red: Variation DBI; Green: Cumulative DBI; Blue: Maximum DBI)
... 20 Figure 4 Sampling sites occupied by Lycosidae, Linyphiidae, Theridiidae, and Oonopidae had higher spider family composition dissimilarity and thus were
divided from other sampling sites into four unique groups through
unweighted pair-group method with arithmetic means (UPGMA) cluster
analysis. Note: figure here only presents families with abundance over 30
viii
individuals. Other 17 Families with abundance less than 30 individuals were
all belonged to the complex community in the cluster analysis. ... 25
Figure 5 Kruskal-Wallis test and pair wised Wilcoxon test indicated that sampling sites dominant by wolf spiders had high cumulative Dynamic Habitat Index (DHI),
Variation DHI, and Variation Dynamic Building/Bareness Index (DBI) while
had low Cumulative DBI. Sampling sites dominant by sheet weaving spiders,
on the other hand, had high cumulative DHI while had low Cumulative DBI,
Variation DHI, and Variation DBI. No statistical significance (Bonferroni
adjusted p-value >0.05) between RS factors at sampling sites dominant by
cobweb spiders and goblin spiders, they both had low cumulative DHI while
had high Cumulative DBI, Variation DHI, and Variation DBI. Upper
boundaries of the rectangles: 25th percentiles; lower boundaries of the
rectangles: 75th percentiles; horizontal bars: median; dotted vertical bars:
upper and lower distribution limits. DHI: Dynamic Habitat Index; DBI:
Dynamic Building/Bareness Index. ... 26
Figure 6 relationship between remote sensing derived environmental factors and sampling site with different spider family composition. DHI: Dynamic
Habitat Index; DBI: Dynamic Building/Bareness Index; Cum: cumulative;
Var: variation. ... 27
ix
Figure 7 Comparison between the proportion of the spider family composition variances explained by remote sensing derived environmental factors (RS
factors) and local vegetation factors (LV factors). (a) RS factors and LV
factors could explain different parts of the spider family composition
variances and together explained 5.91% of total variances. (b) RS factors
could explain more variance (blue circle) than local vegetation factors did
(red circle). SR: species richness; BA: basal area; CV: coverage; DHI:
Dynamic Habitat Index; DBI: Dynamic Building/Bareness Index. ... 28
Figure 8 The proportion of the spider family composition variances explained by remote sensing derived environmental factors at 1, 6, 25, and 100 ha sampling
unit size. At 25 ha spatial scale, remote sensing derived environmental factors
explain the largest percentage of spider family composition variances. ... 29
Figure 9 Model fitness measured by model deviance and pseudo R square of the best- fitted model at 1, 6.25, 25, and 100 ha sampling unit size. The increasing trend
of the model fitness declined at sampling unit size greater than 6.25 ha. ... 33
x
LIST OF TABLES
Table 1 Sites compilation by buffer zones and grids at different spatial scale; CHAO1:
Chao (1984) index estimated species richness; GLM: generalized linear
model; GLMM: generalized linear mix model ... 15
Table 2 Landsat 8 images selected in this study ... 16 Table 3 Pearson correlation coefficient between six remote sensing derived environmental factors; bold characters represent factor pairs with correlation
coefficient over 0.7. Cum_DHI: Cumulative Dynamic Habitat Index;
Min_DHI: Minimum Dynamic Habitat Index; Var_DHI: Variation Dynamic
Habitat Index; Cum_DBI: Cumulative Dynamic Building/Bareness Index;
Max_DBI: Minimum Dynamic Building/Bareness Index; Var_DBI: Variation
Dynamic Building/Bareness Index. ... 18
Table 4 factors included in the analysis and their ecological meanings; EVI: Enhance Vegetation Index; EBBI: Enhanced Built-Up and Bareness Index; DHI:
Dynamic Habitat Indices; DBI: Dynamic Building/Bareness Indices. ... 23
Table 5 Correlation relationship and model AIC between estimated spider species richness and each remote sensing derived environmental factors at 1 ha
sampling unit size. “***”: p-value<0.001; “**”: p-value<0.01; DHI:
xi
Dynamic Habitat Index; DBI: Dynamic Building/Bareness Index. ... 30
Table 6 Correlation relationship between estimated spider species richness and both remote sensing derived environmental factors and local vegetation factors.
“***”: p-value<0.001; SR: species richness; BA: basal area; CV: coverage;
DHI: Dynamic Habitat Index; DBI: Dynamic Building/Bareness Index. ... 31
Table 7 Best-fitted model when modeling estimated spider species richness by remote sensing derived environmental factors at 1 ha, 6.25 ha, 25 ha, and 100 ha
sampling unit size. “***”: p-value<0.001; “**”: p-value<0.01; “*”: p-
value<0.05; “NS”: p-value>0.05; DHI: Dynamic Habitat Index; DBI:
Dynamic Building/Bareness Index. ... 34
1
Introduction
In the recent decades, urban areas expand dramatically (United Nation, 2018), and
urbanization has thought to be one of the leading causes of biodiversity loss (Marzluff,
2001; Li et al., 2005; McKinney, 2008). Therefore, evaluating the influences of
urbanization and the spatial distribution of biodiversity in urban areas are essential keys
to conservation (McKinney, 2008; De Mas et al., 2009; Zhang et al., 2018a).
Spiders are very sensitive to environmental changes both at large and small spatial
scales, which make them an excellent environmental indicator (Pearson, 1994; Fan,
2007; Chapman, 2009; Lowe et al., 2017), and their species richness can be benefitted
by both local and regional scale conservation works (Lowe et al., 2017). Besides, as the
main arthropod predator in the terrestrial ecosystems, spiders can benefit human society
by serving as capable pest controllers (Nyffeler, 1987; Marc et al., 1999; Denno et
al.,2003).
Spider is taxa vary in size and foraging strategy (Uetz et al., 1999; Cardoso et al.,
2011). Therefore, urbanization might cause different effects on different spider families
and functional guilds (Miyashita et al.,1998; Magura et al., 2010; Varet et al., 2014).
Large-sized spiders are often the main victims of urbanization. Natural habitats with
complex vertical vegetation structures where harbor higher prey diversity and biomass
2
are more likely to support high energy cost of large-sized spiders (Miyashita et al.,1998;
Floren and Deleman-Reinhold, 2005; Finch, 2005; Barlow et al., 2007; Hung et al.,
2008; Pinto-Leite et al., 2008; Yekwayo et al., 2016; Setiawan et al., 2016; Lowe et al.,
2017). On the other hands, small-sized space web builders, for example, Theridiidae,
prefer highly urbanized habitats and sometimes become dominant taxon due to their low
energy cost and tolerance to low humidity (Magura et al., 2010). Similar preferences
can also be found in some ground wandering spider taxa (Magura et al., 2010). Orb wed
builder prefer edge habitats or open areas because of their unique spatial needs for orb
web weaving (Miyashita et al.,1998; Petcharad et al., 2016). However, previous studies
also showed that these urban taxa tend to present smaller body size and web width than
those in natural habitats due to the lower prey biomass (Miyashita et al.,1998; Dahirel et
al., 2018). Nevertheless, like most of the wildlife, urbanization could cause negative
influences on overall spider species richness (Miyashita et al.,1998; McKinney, 2008;
Varet et al., 2014).
Remote sensing data can serve as proper tool to monitor the relationship between
broad-scale environmental changes and wildlife species richness (McFeeters, 1996; Liu
et al., 2004; Schowengerdt and Robert, 2007; Mohapatra et al., 2014; Lu et al., 2015).
By satellite or aircraft carried sensors, such as Moderate Resolution Imaging
Spectroradiometer (MODIS), remote sensing can record electromagnetic wave emitted
3
or reflected from the surface instantly (Schowengerdt and Robert, 2007; Lu et al.,
2015). Based on remote sensing data, researchers can derive important environmental
indices to wildlife by calculating the albedo differences between several different
spectral bands reflect from the ground surfaces (McFeeters, 1996; Mohapatra et al.,
2014; Yan et al., 2019).
Vegetation Indices are one of the most commonly used remote sensing derived
environmental indices which are derived from albedo ratio between red (564–580 nm)
and near-infrared (750–950 nm) lights (Liu et al., 2004; Phillips et al., 2010). These
indices are designed to predict the amount of foliage and vegetation density on the
ground (Liu et al., 2004; Phillips et al., 2010) and can also be used on quantifying the
main impacts of urbanization on biodiversity: habitat loss and habitat degradation
(McFeeters, 1996; Li et al., 2005; As-syakur et al., 2012; Mohapatra et al., 2014; Zhang
et al., 2018b).
Based on the albedo ratio, there are several types of formulae to derive vegetation
indices (Liu et al., 2004). Among all vegetation indices, Normalized Difference
Vegetation Index (NDVI) has been the most commonly used vegetation index (Liu et
al., 2004). Derived by red and infrared bands only, NDVI is the most straightforward vegetation index, which can be calculated easily (Liu et al., 2004). However, NDVI also
has several shortages. For example, soil and airborne particle noises can heavily
4
influence the accuracy of NDVI (Liu et al., 2004). Besides, a nonlinear relationship is
presented between the index and real foliage amount. The NDVI saturates at condense
vegetation, such as broadleaf forests (Liu et al., 2004; Phillips et al., 2010; Hobi et al.,
2017). Despite many disadvantages, NDVI is still a sensitive index in habitats with low
to medium vegetation density (Liu et al., 2004; Phillips et al., 2010).
Enhance Vegetation Index (EVI) was then constructed in order to fix the shortages
of NDVI (Hobi et al., 2017; Philpot, 2017). Intaking an additional blue band and two
constants to calibrate soil and practical biases, EVI is more accurate than NDVI for
dense vegetation (Phillips et al., 2010). Therefore, EVI could present linear correlations
to real foliage amount and vegetation density, which make EVI better predictors of
foliage amount than NDVI in tropical and subtropical landscapes (Hobi et al., 2017;
Philpot, 2017). However, no matter how precise a vegetation index is, it does not fully
equal to ground truth foliage amount (Liu et al., 2004; Phillips et al., 2010; Hobi et al.,
2017).
Further, vegetation types and their phenology can be inferred by combining
vegetation indices at different temporal scales (Berry et al., 2007). Dynamic Habitat
Indices (DHIs) are designed to describe vegetation phenology through times base on
vegetation indices (Berry et al., 2007). The DHIs include three different indices,
Cumulative DHI, Minimum DHI, and Variation DHI (Berry et al., 2007). Within a year,
5
Cumulative DHI represents the total amount of foliage of the habitat in the whole year,
Minimum DHI represents the amount of foliage in non-growing season of the habitat,
and Variation DHI is the differences in the amount of foliage between growing and non-
growing season of the habitat (Berry et al., 2007; Andrew et al., 2012; Hobi et al.,
2017). By comparing DHIs in different habitats, vegetation types and their phenology
can be inferred. For example, within one year, temperate deciduous forest presents high
cumulative DHI, high variation DHI, and low minimum DHI because of the dense
foliage in growing season and defoliation in the non-growing season (Hobi et al., 2017).
On the other hand, all three DHIs are low in temperate grassland due to its stable, but
sparse foliage amount (Hobi et al., 2017).
Through vegetation indices and DHIs, researchers can also predict habitat
productivity, which has been thought to be one of the most important environmental
factors to species diversity (Clark et al., 2001; Amthor and Baldocchi, 2001; Liu et al.,
2004; Lu et al., 2015). Habitat productivity is defined as the total inorganic energy fixed
by producers (Amthor and Baldocchi, 2001; Clark et al., 2001). Habitat productivity
hypothesis is an important theory which states that species richness increases with
increasing habitat productivity. However, there are many exceptions which have
rejected habitat productivity hypothesis (McGoff et al., 2013; Vilar et al., 2014; Bicudo
et al., 2016; Zorzal-Almeida et al., 2017). One of the famous cases is eutrophication in
6
aquatic ecosystems, high productivity sometimes leads to explosive growth of bacteria
and algae which consume all the oxygen and cause suffocation to other organisms
(McGoff et al., 2013; Vilar et al., 2014; Bicudo et al., 2016; Zorzal-Almeida et al.,
2017).
Besides species richness, community dissimilarity can also be explained by habitat
productivity (Harrison et al., 2006; Andrew et al., 2012; Johnson and Angeler, 2014). In
habitats with high productivity, community dissimilarity between sampling sites are
larger than those in habitat with low productivity (Harrison et al., 2006; Andrew et al.,
2012; Johnson and Angeler, 2014). Further, Andrew et al. (2012) stated that
dissimilarity is driven by inter-species competition in productive habitats due to the
competition on resources. On the other hand, dissimilarity in less productive habitats are
driven by environmental filtering since only a few strong species are present in harsh
habitats (Andrew et al., 2012).
Although spiders are small-sized carnivores, their species richness and assemblage
composition are heavily influenced by bottom-up effect, prey biomass and vertical
vegetation structure (Greenstone, 1984; Tso, 2003; Fan, 2007; Scharf et al., 2011;
Štokmane and Spuņģis, 2014). Therefore, habitat productivity hypothesis can also apply
on spider species richness (De Mas et al. 2009; Birkhofer and Wolters, 2012; Turrini et
al., 2015). Turrini et al. (2015) even stated that green areas in highly urbanized
7
landscape with high NDVI value can harbor higher spider species richness than that in
rural farmland. Besides species richness, prey diversity to spiders also increases with
increasing productivity at a global scale (Birkhofer and Wolters, 2012). However, in
some cases, spider species richness cannot be well explained or predicted by vegetation
indices (Jiménez-Valverde and Lobo, 2006; Lafage et al., 2014). Lafage et al. (2014)
suggest that the results may be caused by the higher trophic levels that spiders occupied
(Girard et al., 2011). Thus their diversity is not heavily influenced by vegetation
productivity in less productive habitats such as temperate floodplains (Lafage et al.,
2014).
DHIs are good environmental predictors which contain not just spatial, but also
temporal productivity features. To many taxa, habitats with high cumulative, high
minimum productivity, and low productivity variation can harbor higher species
richness, abundance, and larger body size (Coops et al., 2009a; Coops et al., 2009b;
Coops et al., 2009c; Birkhofer and Wolters, 2012; Michaud et al., 2014; Suttidate, 2016;
Hobi et al., 2017; Khlifa et al., 2017; Tiede et al., 2018). However, different taxa and
functional groups might prefer habitats which had different productivity features. To
North American birds, vegetation indices derived DHIs are valid predictors (Hobi et al.,
2017). For forest nesting birds, cumulative DHI has the best predictive power and
presents positive correlation relationship (Hobi et al., 2017). On the other hand,
8
minimum and variation DHI are better predictors for grassland nesting birds, and habitat
had low minimum and high productivity variation harbor higher grassland nesting birds
(Hobi et al., 2017). To spiders, however, the relationship between DHIs and the habitat
preferences of each family or functional guilds was seldom studied.
Habitat preferences of wildlife taxa and functional guilds are not the only reason
which might cause different relationship between habitat productivity and species
richness. Habitat productivity hypothesis is a scale-dependent hypothesis, several
previous studies indicated that the hypothesis is more likely to be applied at broader
spatial scales, while the correlation between habitat productivity and species richness
might become neutral, or even negative and less significant at finer spatial scales
(Coops et al., 2009a; Coops et al., 2009b; Coops et al., 2009c; Birkhofer and Wolters,
2012; Michaud et al., 2014; Suttidate, 2016; Hobi et al., 2017; Khlifa et al., 2017).
Therefore, spatial scale shall also be considered as an important factor when modeling
species richness by habitat productivity. On the other hand, the scale of sampling unit is
also important when using remote sensing data to derived environmental factors (Liu et
al., 2004; De Mas et al. 2009; Phillips et al., 2010). Under remote sensing, differences
within sampling units are shrunk into a single value, which may cause two distinct
habitats had the same environmental factors (Liu et al., 2004; De Mas et al. 2009;
Phillips et al., 2010). To spiders, De Mas et al. (2009) indicated that spider species
9
richness in the mountains in the Mediterranean climate could be well predicted by
NDVI and the area ratio of different vegetation types. They also stated that the proper
spatial scale to predict spider species richness and assemblage composition is in 220-
meter pixel size, and the efficiency decreases at other spatial scales (De Mas et al.
2009).
Not just vegetation indices can quantify the effect of urbanization, building indices,
based on albedos differences between short wave infrared (2000 – 2500 nm) and
thermal infrared (1000 – 1200 nm), are indices sensitive in detecting concrete, asphalt,
and most of the artificial surfaces (McFeeters, 1996; As-syakur et al., 2012; Mohapatra
et al., 2014). Therefore, scientists can map level of urbanization by deriving indices
such as Enhance Build-Up/Bareness Index (EBBI) (As-syakur et al., 2012). However,
only a few studies have focused on the relationship between building indices and
species richness (Shih, 2018). After all, modeling spatial distribution of species richness
at regional scale has become more important to conservation works, especially to those
in urban areas (Turrini et al., 2015; Buchholz et al., 2018; Lowe et al., 2018). Through
burgeoning remote sensing technology, environmental factors which used to cost many
resources to record are now easier to reach (McFeeters, 1996; Liu et al., 2004;
Schowengerdt and Robert, 2007; Mohapatra et al., 2014; Lu et al., 2015). However,
algorithms and conservation strategy were still remaining ambiguous to taxa which
10
attracted less attention, such as spiders (Nyffeler, 1987; Marc et al., 1999; Denno et
al.,2003; McKinney, 2008). Therefore, the present study aimed to test the possibility of
studying species richness and family composition of ground-dwelling arthropods
predator in urbanization landscape by remote sensing data. Further, are spiders in
different families in urbanization landscapes all prefer habitats with the same
productivity features? Also, is the best scale of sampling units for studying spider
diversity in the present study the same to that stated in Lafage et al. (2014) in mountains
in Mediterranean climate (220 m2)? After all, could any conservation strategy be stated
through the present study?
Specifically, the present study aimed to (1) examine differences in spider family
composition in a subtropical urbanization landscape, (2) define habitat preferences of
dominant spider families by remote sensing derived environmental factors (RS factors)
(EVI derived DHIs, and EBBI derived Dynamic Building/Bareness Indices (DBIs)), (3)
explain these spider family composition differences by RS factors, (4) determine
appropriate spatial scale of the sampling units that RS factor(s) can explain the largest
proportion of variance of spider family composition variances, (5) examine
relationships between RS factors and spider species richness in subtropical urbanization
landscape, (6) explore the possibility of modeling spider species richness by RS factors,
and (7) determine best RS factors subset and spatial scale of the sampling units that
11
spider species richness can be modeled by RS factor(s).
12
Methods
GREEnS data set
Global Research & Education on Environment and Society (GREEnS) was a
project supported by the Board of Directors of Tunghai University, and conducted by
Tunghai University from the year 2012 to 2015. One of the main aims of the project was
to test the suitability and validity of studying spider and plants species diversity in
different land use types by using an environmental indicator which has been widely
used in gardening and urban planning, Biotype Area Factor (BAF).
Therefore, a broad urbanization landscape (about 3,500 km2) located in central
Taiwan was selected as the study area of GREEnS project (Figure 1). Within this study
area, a total of 1,232 sampling sites were conducted under 53 patches with four different
land use types, which included 8 school campus, 11 industrial areas, 30 parks, and 4
paddy areas (Appendix 1). Sampling sites was at least 100 m apart and a sampling plot
of 100 m2 quadrat was established for recording both spider and vegetation species
diversity data.
13
Figure 1 The study area located in central Taiwan; each dot represent one sampling site and colors of the dots represent different land use types in Global Research & Education
on Environment and Society project.
Species diversity data
During the growing seasons (March to June) from 2012 to 2015, pitfall traps were
established in the field to record spider species richness and abundance data at every
sampling site and without repeat survey. A set pitfall trap consisted of five 50 ml
centrifuge tubes. Each tube was filled with 30 ml 70% ethanol and left in the field for
14
no more than a week to reach the highest capture efficiency (Schirmel et al., 2010). All
the trapped spider individuals were the removed from the traps and preserved in 70%
ethanol and identified to morphological species level by comparing their genitals. A
total of 1849 spider individuals that belong to 196 morphological species, 27 families
were recorded in the data set (Appendix 2). The present study only included sampling
sites with mature spider individuals and sampling sites with only one juvenile
individual, which was treat as one morphological species at that site. In a total of 1,145
sampling sites were included for further analysis.
On the other hand, vegetation diversity data, including both wood and herb plants,
were surveyed within the 10×10 m sampling plot. For woody plants, individuals with its
diameter at breast height (DBH) larger than 1 cm were identified to species level, and
their DBH were also recorded. For herbs, individuals were identified to species level,
and percentage of coverage of each herb species was also recorded.
In order to determine suitable spatial scale to explain spider species richness and
family composition by RS factors, I rearranged the spider and vegetation diversity data
into four different grid systems, 1ha (100 m diameter buffer of each sampling site), 6.25
ha (250m × 250m grid), 25 ha (500 m × 500 m grid), and 100 ha (1,000 m × 1,000 m
grid) (Table 1). Spider and vegetation diversity data were combined in each grid.
Besides, since no repeat survey was done to both spider and vegetation data.
15
Hence, to correct the possible bias of the data set, I adopted species richness estimation
by using Chao (1984) index (CHAO1), which is a species richness estimator suitable for
estimating invertebrate richness at genus level or lower (Basualdo, 2011). CHAO1 is
derived by the cumulative species richness rarefaction curve and the number of rare
species recorded in sampling sites and is suitable for. The correction method was
applied to spider species richness at every spatial scale (grid size) and to tree and herb
species richness in grids larger than 1 ha (Table 1)
Table 1 Sites compilation by buffer zones and grids at different spatial scale; CHAO1:
Chao (1984) index estimated species richness; GLM: generalized linear model; GLMM:
generalized linear mix model
Unit size
Sample size
CHAO1 Model
Radom factor
N.O. of site(s) per grid
1 ha 1,145
Spider
Richness
GLM None 1
6.25 ha 805
Spider, Wood,
Herb Richness
GLMM
N.O. of sites
per grid
1–9
25 ha 407 1–29
100 ha 186 1–49
16
Remote sensing data
I downloaded five cloud-free images, four from the four seasons of the year 2014
and one from the winter of the year 2013 (Table 2), of the study area (Path: 118, Row:
43) from National Aeronautics and Space Administration (NASA), Landsat 8 mission
level two product online open access (Resolution: 30m × 30m). Based on the images,
Enhance Vegetation Index (EVI) and Enhanced Built-Up and Bareness Index (EBBI)
were calculated through the following formulae (NIR: near infrared band; RED: red
band; BLUE: blue band; SWIR: short wave infrared; TIR: thermal infrared).
EVI = 6 × NIR − RED
NIR + 6 × RED − 7.2 × BLUE + 1 EBBI = SWIR − NIR 10√(SWIR + TIR)
Table 2 Landsat 8 images selected in this study
Date of the image Season represented
December 3rd, 2013 2013 Winter
March 25th, 2014 2014 Spring
August 16th, 2014 2014 Summer
October 19th, 2014 2014 Fall
December 22nd, 2014 2014 Winter
17
Using the EVI and EBBI values in each season, I derived Dynamic Habitat Indices
(DHIs), and Dynamic Building/Bareness Indices (DBIs). Dynamic Habitat Indices
include cumulated annual productivity (Cumulative DHI), minimum annual
productivity (Minimum DHI), and seasonal coefficient (Variation DHI) (Figure 2). On
the other hand, DBIs is consists of cumulated annual build-up/bareness (Cumulative
DBI), maximum annual build-up/bareness (Maximum DBI), and seasonal coefficient
(Variation DBI) (Figure 3).
Both Cumulative DHI and DBI were calculated by integrating EVI and EBBI
value from the winter of the year 2013 to the winter of the year 2014 (Table 3). On the
other hand, Variation DHI and DBI were calculated by the standard deviation of EVI
and EBBI value divided by the average of EVI and EBBI value (Table 3). Minimum
DHI was the lowest EVI value of the year, and Maximum DBI was the largest EBBI
value of the year. All the calculations were done using QGIS version 2.18.15.
However, strong collinearity was found between Cumulative DHI and Minimum
DHI (Pearson correlation coefficient = 0.95), and also between Cumulative DBI and
Maximum DBI (Pearson correlation coefficient = 0.83) (Table 3). Therefore, I excluded
Minimum DHI and Cumulative DBI from the following statistical analyses.
18
Table 3 Pearson correlation coefficient between six remote sensing derived
environmental factors; bold characters represent factor pairs with correlation coefficient
over 0.7. Cum_DHI: Cumulative Dynamic Habitat Index; Min_DHI: Minimum
Dynamic Habitat Index; Var_DHI: Variation Dynamic Habitat Index; Cum_DBI:
Cumulative Dynamic Building/Bareness Index; Max_DBI: Minimum Dynamic
Building/Bareness Index; Var_DBI: Variation Dynamic Building/Bareness Index.
Cum_DHI Min_DHI Var_DHI Cum_DBI Max_DBI Var_DBI
Cum_DHI – – – – – –
Min_DHI 0.95 – – – – –
Var_DHI -0.37 -0.57 – – – –
Cum_DBI -0.21 -0.24 0.20 – – –
Max_DBI -0.14 -0.27 0.37 0.83 – –
Var_DBI -0.33 -0.47 0.52 0.37 0.65 –
19
Figure 2 Dynamic Habitat Indices (DHIs) of the study area derived by Enhance Vegetation Index. (a) cumulative DHI, (b) minimum DHI, (c) variation DHI, (d)
combined DHIs were transformed into RGB color on the map (Red: Variation DHI;
Green: Cumulative DHI; Blue: Minimum DHI)
20
Figure 3 Dynamic Building/Bareness Indices (DBIs) of the study area derived by Enhance Build-up/Bareness Index. (a) cumulative DBI, (b) maximum DBI, (c) variation
DBI, (d) combined DBIs were transformed into RGB color on the map (Red: Variation
DBI; Green: Cumulative DBI; Blue: Maximum DBI)
21
Statistical analyses
All statistical works were done using R version 3.4.3. At all four sizes of sampling
units, eight environmental factors were selected (Table 4), wood species richness, wood
total basal area, herb species richness, and herb coverage were local vegetation factors
(LV factors); and Cumulative DHI, Variation DHI, Cumulative DBI, and Variation DBI
were remote sensing derived environmental factors (RS factors) (Table 4).
Matrix of spider family occupancy at each sampling sites (1 ha scale) was
transformed into Ward’s distance matrix, and went under unweighted pair-group method
with arithmetic means (UPGMA) cluster analysis to group sampling sites based on
spider family composition dissimilarity. Then, Kruskal-Wallis test and pair-wised
Wilcoxon test (Bonferroni adjusted p-value) was adopted to test were RS factors
different between different groups. Together with the results from pair-wised Wilcoxon
test, RS factors preference of dominant spider families in each group and the
significance of studying spider family composition through RS factors were tested by
Bray-Curtis distance based Redundancy Analysis (RDA) and Monte Carlo permutation
test.
Besides, in order to test if RS factors were suitable for explaining spider family
composition differences between sampling sites (1 ha scale), the proportion of
composition variances explained by RS factors and by LV factors through Bray-Curtis
22
distance based RDA were compared. The same comparisons were also made at the 6.25,
25, 100 ha sampling unit size, and the sampling unit size which RS factors could
explain the largest proportion of composition variances was considered as the proper
spatial scale for studying spider family composition with RS factors.
On the other hand, to test the relationship between CHAO1 and RS factors among
sampling sites (1 ha scale), and test did RS factors fit good models to CHAO1, negative
binomial-based generalized linear models (GLM) were conducted between CHAO1 and
each individual RS factor. Also, another two negative binomial-based GLM were
conducted, one with both RS factors and LV factors while the other with RS factors
only. The former was conducted in order to test could any RS factors still reach the 0.05
significance level when modeling CHAO1 with LV factors together. The latter was
conducted for selecting the best-fitted model with highest AIC weight.
Besides, at sampling unit size larger than 1 ha, negative binomial-based
generalized linear mixed models (GLMM) were conducted to model CHAO1 by RS
factors (Table 1). Model AIC and AIC weight of all candidate models with different RS
factors subset combinations were compared to select the best-fitted model for modeling
CHAO1 at each spatial scale. Besides, I set the number of sampling site within a grid as
a random factor to correct its effect on CHAO1. After all, model fitness of the best-
fitted model at each scale was compared by pseudo R square (one minus ration between
23
null deviance and residual deviance) between the predicted value generated from the
model and CHAO1. The spatial scale that best-fitted model had the highest pseudo R
square was considered the proper spatial scale of modeling CHAO1 by RS factors.
Table 4 factors included in the analysis and their ecological meanings; EVI: Enhance Vegetation Index; EBBI: Enhanced Built-Up and Bareness Index; DHI: Dynamic
Habitat Indices; DBI: Dynamic Building/Bareness Indices.
Type Factors Calculation Ecological meaning
Local vegetation
Wood species richness — Horizontal vegetation structure
Wood basal area π × (DBH 2 )
2
144
Vegetation abundance
Herb species richness — Horizontal vegetation structure
Herb coverage — Vegetation abundance
Remote sensing
Cumulative DHI ∫ EVI
2013 Dec 2014 Dec
Annual total productivity
Variation DHI ∂(EVI)
μ(EVI) Annual productivity variation
Cumulative DBI ∫ EBBI
2013 Dec 2014 Dec
Annual total Building/Bareness
Variation DBI ∂(EBBI) μ(EBBI)
Annual Building/Bareness
Variation
24
Results
Spider family composition differences in the subtropical urbanization landscape
were mainly determined by the occupancy of four dominant spider families, which were
wolf spiders (Lycosidae), sheet weaving spiders (Linyphiidae), goblin spiders
(Oonopidae), and cobweb spiders (Theridiidae) (Figure 4). Based on the proportion of
the occupancy of the four dominant spider families, sampling sites were divided into
five groups: wolf spiders dominance (occupancy>0.5), sheet weaving spiders
dominance, goblin spiders dominance, cobweb spiders dominance, and complex
community which included other 23 spider families with low species richness and
abundance (Figure 4).
Remote sensing derived environmental factors (RS factors) at sampling sites in
different groups were significantly different (Figure 5). Thus, habitat preferences of
these dominance spider families could be stated by the RS factors (Figure 6). Wolf
spiders preferred habitats which had high but unstable annual productivity and low but
also unstable annual building/bareness (Figure 6). Sheet weaving spiders, on the other
hand, preferred habitats which had high and stable annual productivity and low and also
stable annual building/bareness (Figure 6). Goblin spiders and cobweb spiders,
however, preferred similar habitats which had low and unstable annual productivity
25
while had high and unstable annual building/bareness (Figure 6). Beside the four
dominant families, most of the less-abundant spider families were recorded in habitats
which had low but stable annual productivity and medium but stable annual
building/bareness.
Figure 4 Sampling sites occupied by Lycosidae, Linyphiidae, Theridiidae, and Oonopidae had higher spider family composition dissimilarity and thus were divided
from other sampling sites into four unique groups through unweighted pair-group
method with arithmetic means (UPGMA) cluster analysis. Note: figure here only
presents families with abundance over 30 individuals. Other 17 Families with
abundance less than 30 individuals were all belonged to the complex community in the
26
cluster analysis.
Figure 5 Kruskal-Wallis test and pair wised Wilcoxon test indicated that sampling sites dominant by wolf spiders had high cumulative Dynamic Habitat Index (DHI), Variation
DHI, and Variation Dynamic Building/Bareness Index (DBI) while had low Cumulative
DBI. Sampling sites dominant by sheet weaving spiders, on the other hand, had high
cumulative DHI while had low Cumulative DBI, Variation DHI, and Variation DBI. No
statistical significance (Bonferroni adjusted p-value >0.05) between RS factors at
sampling sites dominant by cobweb spiders and goblin spiders, they both had low
cumulative DHI while had high Cumulative DBI, Variation DHI, and Variation DBI.
27
Upper boundaries of the rectangles: 25th percentiles; lower boundaries of the rectangles:
75th percentiles; horizontal bars: median; dotted vertical bars: upper and lower
distribution limits. DHI: Dynamic Habitat Index; DBI: Dynamic Building/Bareness
Index.
Figure 6 relationship between remote sensing derived environmental factors and sampling site with different spider family composition. DHI: Dynamic Habitat Index;
DBI: Dynamic Building/Bareness Index; Cum: cumulative; Var: variation.
RS factors could serve as significant constrained factors in explaining spider
family composition variances and could explain variances which LV factors could not
28
explain (Figure 7). Further, among the variances explained by both types of factors
(5.91% of the total variances), larger proportion of variances were explained by RS
factors than LV factors did (RS factors: 41%; LV factors: 35%; both: 24%) (Figure 7).
Also, the proportion of spider family composition variances explained by RS factors
increased about 10% from 1 ha spatial scale to 25 ha spatial scale while declining at the
spatial scale larger than 25 ha spatial scale (25 ha: 53%; 100 ha: 50%) (Figure 8).
Therefore, 25 ha spatial scale was the proper scale for using RS factors on explaining
spider family composition.
Figure 7 Comparison between the proportion of the spider family composition variances explained by remote sensing derived environmental factors (RS factors) and
local vegetation factors (LV factors). (a) RS factors and LV factors could explain
different parts of the spider family composition variances and together explained 5.91%
29
of total variances. (b) RS factors could explain more variance (blue circle) than local
vegetation factors did (red circle). SR: species richness; BA: basal area; CV: coverage;
DHI: Dynamic Habitat Index; DBI: Dynamic Building/Bareness Index.
Figure 8 The proportion of the spider family composition variances explained by remote sensing derived environmental factors at 1, 6, 25, and 100 ha sampling unit size.
At 25 ha spatial scale, remote sensing derived environmental factors explain the largest
percentage of spider family composition variances.
30
Estimated spider species richness (CHAO1) was higher in habitats which had high
annual productivity while lower in habitat with high productivity variation, high annual
building/bareness, and high building/bareness variation and all four RS factors were
significantly correlated to CHAO1 (Table 5). Besides, although the model with local
vegetation factors (LV factors) added fitted CHAO1 better (AIC with LV factor: 3887.3;
AIC without LV factor: 3982.7), RS factors could still found significantly correlated to
CHAO1, especially Cumulative DHI (Table 6).
Table 5 Correlation relationship and model AIC between estimated spider species richness and each remote sensing derived environmental factors at 1 ha sampling unit
size. “***”: p-value<0.001; “**”: p-value<0.01; DHI: Dynamic Habitat Index; DBI:
Dynamic Building/Bareness Index.
Factor Trend Significance AIC
Cumulative DHI + *** 3985
Variation DHI – *** 4129
Cumulative DBI – ** 4166
Variation DBI – *** 4133
31
Table 6 Correlation relationship between estimated spider species richness and both remote sensing derived environmental factors and local vegetation factors. “***”: p-
value<0.001; SR: species richness; BA: basal area; CV: coverage; DHI: Dynamic
Habitat Index; DBI: Dynamic Building/Bareness Index.
Factor type Factor Trend Significance
Local vegetation
Wood SR + ***
Wood BA +
Herb SR + ***
Herb CV + ***
Remote sensing
Cumulative DHI + ***
Variation DHI –
Cumulative DBI +
Variation DBI –
AIC 3887.3
RS factors subsets significantly correlated to CHAO1 at 1, 6.25, and 25 ha
sampling unit size (Table 7). The best-fitted model at 1 ha spatial scale consisted of
Cumulative DHI and Variation DBI (Weight AIC: 0.383), the former was positively
32
correlated to CHAO1 while the latter was negatively correlated to CHAO1 (Table 7). At
6.25 ha spatial scale, Cumulative DHI was the significant factor which positively
correlated to CHO1 and fitted CHAO1 better than other candidate models did at the
same spatial scale (AIC weight: 0.267) (Table 7). On the other hand, The best-fitted
model at 25 ha spatial scale consisted of Cumulative DHI and Variation DHI (Weight
AIC: 0.339), the former was positively correlated to CHAO1 while the latter was
negatively correlated to CHAO1 (Table 7). At 100 ha spatial scale, the model with
Cumulative DHI alone fitted CHAO1 better than other candidate models did at the same
scale. However, Cumulative DHI did not reach the 0.05 significance level (Table 7).
Further, model fitness of the best-fitted model at each spatial scale was different, it
increased from 1 ha to 100 ha spatial scale (Pseudo R square at 1 ha: 0.103; 6.25 ha:
0.506; 25 ha: 0.595; 100 ha: 0.622) (Figure 9). However, the increasing trends of the
model fitness declined at sampling unit size greater than 6.25 ha (Figure 9). After all,
considered both model fitness and factor significance of the best-fitted model at each
spatial scale, sampling unit size between 6.25 and 25 ha were the proper sampling unit
size for modeling CHAO1 by RS factors, and Cumulative DHI could serve as an
important RS factor.
33
Figure 9 Model fitness measured by model deviance and pseudo R square of the best- fitted model at 1, 6.25, 25, and 100 ha sampling unit size. The increasing trend of the
model fitness declined at sampling unit size greater than 6.25 ha.
34
Table 7 Best-fitted model when modeling estimated spider species richness by remote sensing derived environmental factors at 1 ha, 6.25 ha, 25 ha, and 100 ha sampling unit size. “***”: p-value<0.001; “**”: p-value<0.01; “*”: p-value<0.05; “NS”: p-value>0.05; DHI: Dynamic Habitat
Index; DBI: Dynamic Building/Bareness Index.
Factor
1 ha 6.25 ha 25 ha 100 ha
Trends Significance Trends Significance Trends Significance Trends Significance
Cumulative DHI + *** + *** + ** + NS
Variation DHI - *
Cumulative DBI
Variation DBI – **
Weight(AIC) 0.383 0.267 0.339 0.143
35
Discussion
Results of this study suggest spider communities in the subtropical urbanization
landscape were mainly dominated by wolf spiders, sheet weaving spiders, cobweb
spiders, and goblin spiders. RS factors could significantly explain the composition
variances and define the habitat preferences of each dominant spider family. Besides,
RS factors could explain the largest proportion of spider family composition variances
at the 25 ha spatial scale. On the other hand, RS factors were found significantly
correlated to CHAO1, and only Cumulative DHI was positively correlated to CHAO1.
Also, the model with Cumulative DHI and Variation DHI fitted CHAO1 the best at the
spatial scales between 6.25 to 25 ha.
Relationship between animal composition assemblage and nature resources is
thought to be driven by two mechanisms, inter-species competition, and environmental
filtering (Andrew et al., 2012; Nakadai et al., 2018). In habitats which are lack of
resources, the effect of environmental filtering is usually larger than competition and
only a few species/taxa of animal can survive become dominant (Andrew et al., 2012).
Therefore, the present study speculates that the differences in ground-dwelling spider
composition in the present study may also be driven by environmental filtering due to
low and unstable productivity in urban areas and its surrounding landscapes. Similar
36
result has been also reported in other invertebrates like butterflies (Andrew et al., 2012).
The main differences which cause those dominant spider families preferred
different habitats shall be their foraging behavior and morphological traits. Most of the
sheet weaving spiders are web builders in small to medium size, which weaves 2-
dimensional sheet web and can hunt for their prey efficiently and control most of the
agricultural pests (Uetz et al., 1999; Cardoso et al., 2011; Yin et al., 2011). However,
web building is also a foraging strategy which costs a lot of energy (Uetz et al., 1999;
Cardoso et al., 2011). Also, sheet weaving spiders need to anchor their web to
vegetation (Cardoso et al., 2011). Therefore, based on the special needs and high energy
cost of web building, sheet weaving spiders need habitats with high and stable year-
round productivity. Therefore, sheet weaving spiders can serve as a good environmental
indicator on monitoring urban greenness due to their habitat preference.
Wolf spiders are ground wandering spiders in medium to large size (Uetz et al.,
1999; Cardoso et al., 2011), they prefer habitat with higher productivity to support the
high energy cost caused by their larger body. However, rather than hunting by web in a
specific spot, wolf spiders have good mobility, which makes them able to travel for a
distance for hunting (Uetz et al., 1999; Cardoso et al., 2011). Their mobility also gives
them flexibility when facing disturbance (Yu et al., 2002). which may make them
possible to dwell in habitats with unstable productivity, such as paddy fields and
37
wetlands.
Goblin spiders are small-sized ground wandering spiders smaller than wolf spiders
do (Cardoso et al., 2011). Therefore, the energy cost of goblin spiders are low. Further,
the exoskeleton in some of the goblin spiders is thickened and skeletonized (Simon,
1890). Altogether, goblin spiders can endure small and dry habitats in the highly
urbanized area such as small vegetated patches or street trees on asphalt and concrete.
Like sheet weaving spiders, most of the cobweb spiders are also web builders in
small to medium size, while they weave a 3-dimensional irregular web for catching
preys (Uetz et al., 1999; Cardoso et al., 2011). However, cobweb spiders have been
reported as edge or urban species in several cases (Miyashita et al.,1998; Magura et al.,
2010). They can anchor their web on artificial structures and forage on the sufficient
insects lured by human activities such as street lights (Magura et al., 2010).
Overall, although spiders are small invertebrate predators, the correlation
relationships between CHAO1 and DHIs are the same to those in vertebrates and
herbivorous invertebrates, such as birds (Coops et al., 2009a; Coops et al., 2009b;
Coops et al., 2009c; Hobi et al., 2017) and butterflies (Andrew et al., 2012). Therefore,
ground-dwelling spiders in subtropical urbanization landscape follow habitat
productivity hypotheses. However, the hypothesis was not supported in spiders in
temperate floodplains landscapes (Lafage et al., 2014) although two landscape types
38
share similar EVI value (maximum about 0.5). Possible reason for this finding is these
two studies were conducted at different spatial scales, size of the study area in Lafage et
al. (2014) was 7.4 ha while that is the present study could up to more than 1000 ha.
Habitat productivity hypothesis is more likely to be agreed at regional scale (Coops et
al., 2009a; Coops et al., 2009b; Coops et al., 2009c; Birkhofer and Wolters, 2012;
Michaud et al., 2014; Suttidate, 2016; Hobi et al., 2017; Khlifa et al., 2017). Besides,
the RS factors in the present study not just contain vegetation indices, but also building
indices. However, like the previous study (Shih, 2017), building indices cannot serve as
factors which are as good as vegetation indices when modeling species richness.
RS factors could explain more variance to both CHAO1 and spider family
composition, and agree with the results reported in previous studies, which stated that
productivity related factors are more significant at regional scale (Coops et al., 2009a;
Coops et al., 2009b; Coops et al., 2009c; Birkhofer and Wolters, 2012; Michaud et al.,
2014; Suttidate, 2016; Hobi et al., 2017; Khlifa et al., 2017). However, the best spatial
scale, which RS factors could explain the largest proportion of both CHAO1 and spider
family composition variances, in the present study was larger than 4.84 ha (220 m2)
suggested by De Mas et al. (2009). Main reason for this difference may be the area ratio
of different land covers in the study area. Study area of De Mas et al. (2009) was
located in mountains (altitude range: 800 m–2500 m), which were landscapes with
39
complex land covers with various vegetation types at different altitudes (Lasseur et al.,
2018). On the other hand, land covers in the urbanization landscape in the present study
were simpler than those in mountains, and land cover of most of the sampling sites in
the present study was large-sized paddies and urban areas. Therefore, large-sized
sampling grids in landscapes with complex land covers may mix multiple land covers
with different spider composition and obscure the relationships of RS factors with
spider species richness and family composition. Overall, the present study suggests that
sampling unit with its size between 6.25 to 25 ha is the best when studying both spider
species and family composition in subtropical urbanization landscape.
Although DHIs was designed to describe vegetation phenology (Berry et al.,
2007), the present study found they can also be used in describing urbanized land cover
types. Parks and school (green area) are land cover which had high cumulative DHI,
high minimum DHI, and low variation DHI; paddies, on the other hand, had high
cumulative DHI, low minimum DHI, and high variation DHI while urbanized area had
low cumulative, minimum, and variation DHI. However, land covers with high
minimum DHI, but low cumulative and variation DHI, which usually represent
temperate grass land in previous studies (Coops et al., 2009a; Coops et al., 2009b;
Coops et al., 2009c; Andrew et al., 2012), cannot be found in the study areas and thus
cause strong collinearity between cumulative and minimum DHI. Therefore, the present
40
study suggests that when studying species diversity in tropical or subtropical landscapes
minimum DHI can be skipped, and represented by cumulative DHI, which contains
more temporal information of primary productivity.
The spider distribution data set in the present study was previously used in a study
(Huang et al., 2015) examing relationship between spider species richness and other
environmental index, Biotope Area Factor (BAF), which is defined as ratios of
ecologically effective land covers (e.g. green areas) among landscapes (Butchar et al.,
2010). BAF has been a commonly used index in urban planning in Europe, and is
derived by scoring every surface material at sampling sites mainly based on their water
permeability (e.g. 5 points for soil while 1 point for asphalt) (Huang et al., 2015).
However, BAF do not present any significant relationship to neither spider, nor plant
species richness (Huang et al., 2015). Therefore, the present study suggests that indices
of habitat productivity shall also be considered when using BAF or other similar
indices.
Previous studies indicated that proper managements in urban green areas, such as
parks, are critical to spider conservation (Buchholz et al., 2018; Lowe et al., 2018).
However, the present study showed that not just green areas, but also small vegetated
patches and even street trees can serve as important habitats to many different spider
families. In Taiwan, these small habitats have often undergone pesticide and herbicide
41
usage in order to keep the urban clean and green. Although the target of pesticide and
herbicide are not spiders, the toxic chemicals can still be consumed by spiders through
bioconcentration effect (Girard et al., 2011; Gill & Garg, 2014; Maurya & Malik, 2016)
and lead to their death by harming their nerve systems (Samu &Vollrath, 1992).
Therefore, beside the effect of low productivity, the use of toxic chemical may also be
one of the main cause of low spider species richness and abundance in highly urbanized
area.
The present study offers methods of modeling the spatial distribution of small
carnivorous invertebrate species richness in the urbanization landscape by a very large
scale monitoring scheme and remote sensing data. In conservation works, the present
study suggests that maintaining high and stable habitat productivity, such as parks or
school campuses, can help maintain high spider species richness. On the other hand,
reduce the amount of pesticide and herbicide usage on street trees and small vegetated
patches in highly urbanized areas, can also help on increasing urban spider species
richness and abundance.
42
References
Amthor, J. S., & Baldocchi, D. D. (2001). Terrestrial higher plant respiration and net
primary production. In: Roy, J., Saugier, B., Mooney, A. H., (eds) Terrestrial
Global Productivity, 33–59, Academic Press, New York.
Andrew, M. E., Wulder, M. A., Coops, N. C., & Baillargeon, G. (2012). Beta‐diversity
gradients of butterflies along productivity axes. Global Ecology and Biogeography,
21(3), 352–364.
As-syakur, A., Adnyana, I., Arthana, I. W., & Nuarsa, I. W. (2012). Enhanced built-up
and bareness index (EBBI) for mapping built-up and bare land in an urban area.
Remote Sensing, 4(10), 2957–2970.
Barlow, J., Gardner, T. A., Araujo, I. S., Ávila-Pires, T. C., Bonaldo, A. B., Costa, J. E.,
& Hoogmoed, M. S. (2007). Quantifying the biodiversity value of tropical primary,
secondary, and plantation forests. Proceedings of the National Academy of
Sciences, 104(47), 18555–18560.
Basualdo, C. V. (2011). Choosing the best non-parametric richness estimator for benthic
macroinvertebrates databases. Revista de la Sociedad Entomológica Argentina,
70(1–2), 27-38.
Berry, S., Mackey, B., & Brown, T. (2007). Potential applications of remotely sensed
43
vegetation greenness to habitat analysis and the conservation of dispersive fauna.
Pacific Conservation Biology, 13(2), 120–127.
Bicudo, D. C., Tremarin, P. I., Almeida, P. D., Zorzal-Almeida, S., Wengrat, S.,
Faustino, S. B., & Morales, E. A. (2016). Ecology and distribution of Aulacoseira
species (Bacillariophyta) in tropical reservoirs from Brazil. Diatom Research,
31(3), 199–215.
Birkhofer, K., & Wolters, V. (2012). The global relationship between climate, net
primary production and the diet of spiders. Global Ecology and Biogeography,
21(2), 100–108.
Buchholz, S., Hannig, K., Möller, M., & Schirmel, J. (2018). Reducing management
intensity and isolation as promising tools to enhance ground-dwelling arthropod
diversity in urban grasslands. Urban Ecosystems, 21(6), 1139–1149.
Cardoso, P., Pekár, S., Jocqué, R., & Coddington, J. A. (2011). Global patterns of guild
composition and functional diversity of spiders. PloS one, 6(6), e21710.
Chapman, A. D. (2009). Numbers of living species in Australia and the world.
Canberra, Australia: Australian Government.
Clark, D. A., Brown, S., Kicklighter, D. W., Chambers, J. Q., Thomlinson, J. R., & Ni, J.
(2001). Measuring net primary production in forests: concepts and field methods.
Ecological Applications, 11(2), 356–370.
