국제 교역의 증가는 다양한 외래 곤충(Exotic insect)들의 국내 유입을 촉진해왔다 (Hong et al., 2012). 또한 최근 빠르게 변화하고 있는 기후는 유입된 외래 곤충들이 새로운 지역에서 정착할 수 있는 가능성을 증가시켰다(Hellmann et al., 2008). 유입된 외래 곤충은 종종 정착에 성공한 후 확산되어 생태적, 경제적으로 부정적인 영향을 야기하는 침입 곤충(Invasive insect) 으로 해충화되어 국내에서 문제를 일으키고 있다(Baek et al., 2019;Kim et al., 2019;Namgung et al., 2020). 이에 환경부에서는 미국선녀벌레, 갈색날개매미충, 꽃매미, 등검은말벌 등의 침입 곤충을 생태계교란생물로 지정하고 모니터링과 관리를 위한 연구를 수행하고 있다(https://kias.nie.re.kr/home/main/main.do).

2023년 5월 22일~23에 농림축산검역본부, 환경부, 문화재청, 산림청 등 유관기관의 협업으로 서울시 강남구 주택에서 외래 흰개미에 대한 합동 역학조사를 실시하였다(https://www.cha.go.kr/newsBbz/selectNewsBbzView.do?newsItemId=155 704127&sectionId=b_sec_1&mn=NS_01_02). 조사 결과 실내 문틀에서 여왕을 포함한 총 159마리의 개체가 확인되었다. 이 종은 Cryptotermes domesticus Haviland (Blattodea: Kalotermitidae)로 확인되었으며(이하 (가칭)마른나무흰개미), 감염된 목재나 가구를 통해 유입된 후 실내에서 군체를 형성한 것으로 보인다(https://www.cha.go.kr/). 국내에서 이번에 발견된 마른 나무흰개미의 원산지는 인도차이나 반도와 인도네시아로 추정되며, 이 밖에도 중국 남부, 대만, 일본, 호주 북부, 폴리네시아 제도, 파나마 등 열대･아열대 지역에서 발견된 기록이 있다 (Scheffrahn and Krecek, 1999;Evans et al., 2013). 침입 흰개미류로 분류되는 종들은 목재를 먹이로 삼는 동시에 이를 둥지로 삼으며, 번식 능력이 뛰어나다(Evans, 2010). 내부에 군체를 형성하는 특징 때문에 Kalotermitidae과의 많은 흰개미 종들이 목재 및 목재 가공품들을 통해 유입될 수 있으며, 이에 따른 유입압(Propagule pressure)이 상대적으로 높은 것으로 평가된다 (Evans, 2010). 2021년에 국내에서 발견된 외래 흰개미도 Kalotermitidae과의 Glyptotermes nakajimai Morimoto로 전라남도 완도군 여서도에서 발견된 바 있다(Shim et al., 2021).

마른나무흰개미는 살아있는 나무를 비롯하여 목재, 가구, 목조건물 등에서 발견된 기록이 있어 국내에 정착 후 확산시에는 잠재적인 해충이 될 가능성이 있다(Himmi et al., 2021). 중국에서는 하이난성에 마른나무흰개미가 침입하였으며, 목재를 가해하는 심각한 해충으로 보고하였다(Huang et al., 2009). 따라서 국내에서도 마른나무흰개미의 확산을 조기에 방지하기 위해 대대적인 발생 현황 조사가 필수적이다. 종 분포 모델 (Species distribution model; SDM)은 발생 기록을 바탕으로 대상 종의 생태적 지위를 추정하는 방법으로써, 침입종의 정착 가능성을 평가하는 데에 널리 이용되고 있다(Elith and Leathwick, 2009). 종 분포 모델로 예측된 서식 적합성 지도는 조사 침입종의 향후 관리 방향에 대해 결정할 수 있는 과학적 근거를 마련 해 줄 수 있을 것이다. 이에 본 연구에서는 마른나무흰개미의 전세계 발생 기록을 이용하여 이들이 정착하기에 적합한 기후 조건을 탐색하는 것을 목표로 하였다. 또한 국내에서 마른나무 흰개미의 기후적합성을 평가하고, 정착 위험이 높은 지역을 파악하고자 하였다.

재료 및 방법

분포자료

마른나무흰개미의 기후 적합성을 평가하기 위해 전세계 분포 자료를 수집하였다. Global Biodiversity information Facility (GBIF; https://www.gbif.org/)에 등재된 134개의 기록을 수집 하였고, 이 중 정확한 지리좌표가 제공되지 않거나 중복되는 지점들을 제외한 35개의 분포 자료를 정리하였다. 또한 문헌 자료로부터 12개 지점의 분포 자료를 추가하여 총 47개 지점의 경위도 좌표로 구성된 분포 데이터를 구축하였다(Indrayani et al., 2017;Himmi et al., 2021)(Fig. 1). 생태 지위 추정을 위해 의사 부재자료를 생성하였다(Barbet-Massin et al., 2012). 의사부재 자료는 대상종의 존재 지점에 대비하여 잠재적으로 부재한 지점(혹은 평균적인 조건을 가지는 지점)의 환경 요인에 대한 차이를 모델링 하는데 이용된다. 잠재적으로 기후적합성이 높을 것으로 예상되는 지점에서의 임의부재자료 추출을 방지하기 위해 47개의 존재 자료로 생성한 볼록다각형(Convex hull)을 이용하였다. 볼록다각형으로부터 2˚(2-decimal-degree)의 버퍼를 생성한 후에 버퍼 바깥쪽 지역에서 1000개의 의사부재자료를 임의로 추출하였다. 분포 자료 생성 과정을 10번 반복하여 47개 존재 자료와 1000개의 의사부재자료가 포함된 총 10개의 분포데이터를 생성하였다. 분포데이터 처리 과정은 R 4.2.2의 sp 패키지를 이용하여 수행하였다(R Core Team, 2021).

생물기후변수

마른나무흰개미의 기후적합성을 추정하기 위해 WorldClim 에서 제공하는 19개의 생물기후변수(Bioclimatic variable)(10 min resolution)를 이용하였다(https://www.worldclim.org/). 입력변수의 다중공선성(Multicollinearity)을 피하고 해석의 용이성을 위해 변수 간의 상관 관계를 검토하였다. 변수 간의 상관지수(Pearson’s r)가 0.85이상인 경우 두 변수 중 1개의 변수만 선택하고 나머지는 분석에서 제외하였다. 이 과정을 통해서 최종적으로 총 10개의 생물기후변수를 모델 구축에 활용하였다. 모델에 입력한 변수는 다음과 같다: 평균 일교차(Bio2), 온도 등온성(Bio3), 연평균 온도 범위(Bio7), 가장 습한 분기의 평균 온도(Bio8), 가장 건조한 분기의 평균 온도(Bio9), 가장 습한 달의 강수량(Bio13), 가장 건조한 달의 강수량(Bio14), 강수량 계절 변동성(변동 계수, Coefficient of Variation) (Bio15), 가장 따듯한 분기의 강수량(Bio18), 가장 추운 분기의 강수량(Bio19).

기후 적합성 모델링

분포 자료의 지리 좌표를 이용하여 각 지점(존재, 의사부재)에 해당하는 10개의 생물기후변수 값을 추출하였다. 이 변수들을 입력 변수로 하고 존재/의사부재 여부를 반응변수로 하여 종 분포 모델을 구동하였다. 존재/의사부재는 1/0값으로 입력하였으며, 모델 알고리즘으로는 Flexible Discriminant Analysis (FDA), Multivariate Adaptive Regression Splines (MARS), Gradient Boosting Machine (GBM), Random Forest (RF)을 이용하였다(Thuiller et al., 2009). 따라서 10개의 데이터세트와 4 개의 알고리즘을 이용하여 총 40개의 모델을 구동하였다. 각 데이터 세트의 70%를 모델 훈련에, 30%를 모델 평가에 이용하였 다. 모델 구동은 R 4.2.2의 Biomod2 패키지에서 이루어졌으며, 사용자가 모델 구동 시 조정해야 하는 하이퍼파라미터 값들은 패키지에서 제공하는 기본값을 사용하였다. 각 모델들은 테스트 데이터를 이용하여 Area Under the Receiver Operating Characteristic Curve (AUC)와 True Skill Statistics (TSS)로 평가되었다(Hao et al., 2019). AUC는 모델이 모든 평가데이터로부터 존재 지점을 얼마나 잘 구별하는지를 측정하는 지표로 0~1의 값을 가지며, 0.5와 1은 각각 무작위 예측과 완벽한 예측을 의미한다. TSS는 존재 지점과 의사부재 지점을 모두 정확하게 예측하는지를 측정하는 지표로 혼동행렬에서 산출되는 민감도(Sensitivity)와 특이도(Specificity)의 합에서 1을 뺀 값으로 산출된다. 따라서 1은 완벽한 모델, 0은 무작위 예측을 의미 한다. 우수한 성능(AUC > 0.8, TSS > 0.7)을 보인 모델들의 결과(존재 확률)를 평균하여 앙상블 예측 모델을 개발하였다 (Thuiller et al., 2009). 최종 앙상블 모델에서 생물기후변수들의 상대적 기여도는 순열 중요도(Permutation importance) 평가를 10회 반복하여 산출되었다. 최종 앙상블 모델을 이용하여 마른나무흰개미의 정착 가능성을 산출하고 지리적 공간에 투영하였다. 국내에서의 기후 적합성은 30 sec의 공간해상도를 생물기후변수로 예측하고 투영하였다.

결과 및 고찰

마른나무흰개미의 발견 기록들을 정리한 결과 대부분 동남 아시아와 태평양의 섬들에 분포하는 것으로 확인되었다(Fig. 1). 대부분의 분포는 섬이나 연안에 집중되는 특징을 보여 침입 후 내륙으로의 자연적 확산이 어렵거나 서식이 적합하지 않을 수 있을 것으로 추정된다.

서식지 적합성 예측 모델에서 AUC (mean ± S.D.)는 RF 모델이 가장 우수하였다: FDA (0.987 ± 0.0047); GBM (0.985 ± 0.0186); MARS (0.944 ± 0.0352); RF (0.989 ± 0.0143). 반면에 TSS (mean ± S.D.)는 FDA 모델이 사용한 알고리즘 중에서 가장 우수하였다: FDA (0.890 ± 0.0432); GBM (0.839 ± 0.0709); MARS (0.850 ± 0.0786); RF (0.847 ± 0.0424). 하지만, 40개의 개별 모델들은 모두 우수한 예측 성능을 보였기 때 문에 모든 모델 결과 값들에 가중치를 부여하지 않고 평균하여 앙상블 모델을 개발하였다. 최종 앙상블 모델에서 연평균온도 범위(Bio7)가 모델 성능에 가장 큰 영향을 미치는 것으로 나타났다(Fig. 2). 따라서 마른나무흰개미는 연평균기온의 편차가 상대적으로 작은 해양성 기후에서 더 잘 정착할 수 있는 것으로 추정되었다(Fig. 3). 이는 현재까지 보고된 분포가 해안 주변으로 한정된다는 결과가 반영된 것으로 보인다. 그 다음으로는 연 교차와 일교차가 비슷한 수준으로 유지되는 지역에서의 기후 적합성이 높았다(Fig. 2). 따라서 사계절이 뚜렷한 온대 지방에서의 기후 적합성이 상대적으로 낮을 것으로 추정되었다.

최종적으로 지리적 공간에 투영된 기후적합성 예측 지도는 동남아시아와 태평양 도서 지역이 마른나무흰개미의 정착에 적합하다는 것을 보여주었다(Fig. 4). 본 연구에서는 정확한 지리 좌표의 부재로 파나마와 트리니다드 토바고에서의 발견 기록을 모델에 입력하지 않았다. 하지만, 최종 모델의 결과는 남미 연안 지역과 카리브해 지역에서의 기후적합성이 상대적으로 높을 것으로 예측하였고, 따라서 이는 실제 발생과 일치하는 결과라고 할 수 있다(Araujo, 1970;Evans et al., 2013).

국내의 마른나무흰개미 기후적합성은 전반적으로 낮을 것으로 예상되었다. 국내로만 한정 했을 때는 제주 지역이 다른 지역에 비해 기후적합성이 약간 높았으나, 원산지에 비해는 현저히 낮은 수준인 것으로 보인다(Fig. 4). 최근 국내에서 발견된 마른나무흰개미는 연안에서 멀리 떨어진 서울 강남에서만 발생한 것으로, 항구를 통해 침입한 개체군이 정착한 후 서울지역으로 확산한 것으로 보기는 어려울 것으로 판단된다. 따라서 이번에 발견된 군체는 목재를 통해 주거지로 유입된 개체들이 비교적 환경 조건의 편차가 작은 실내에서 번식에 성공하여 형성 된 것으로 추정된다. 하지만 마른나무흰개미의 생활사가 명확 하게 밝혀지지 않았고, 분포 기록이 한정되어 있는 만큼 이 종에 대한 모니터링이 항구지역을 중심으로 지속해서 이루어져 야 할 것으로 보인다. 또한 최근의 기후변화에 따라 국내 기후의 아열대화가 진행되고 있는 만큼 국내에서 마른나무흰개미의 정착과 확산 가능성을 염두한 정밀 역학 조사가 필요하다.