1. 서 론

뒤영벌(genus Bombus)은 온대 및 한대 생태계의 핵심 화분매개곤충으로, 야생식물의 수분과 농작물 생산에 중요한 역할을 담당한다(Cameron and Sadd 2020). 특히 뒤영벌은 저온 환경에서도 활동할 수 있는 내한성(Heinrich 1979)과 진동수분(buzz pollination) 능력(De Luca and Vallejo-Marín 2013)을 지녀, 다른 화분매개자가 활동하기 어려운 고위도·고산 환경에서 대체 불가능한 생태적 기능을 수행한다(Goulson 2010). 그러나 최근 전 세계적으로 뒤영벌의 개체군 감소와 분포역 축소가 보고되고 있으며, 그 주요 원인으로 기후변화, 서식지 파괴, 농약 사용, 병원체 확산 등이 지목되고 있다(Goulson et al. 2008;Cameron and Sadd 2020).

기후변화는 뒤영벌 쇠퇴의 특히 중요한 요인으로 주목받고 있다. 북아메리카와 유럽의 대륙규모 분석에 따르면, 뒤영벌은 기후 온난화에 대응하여 분포역의 남쪽 경계가 북상하는 반면 북쪽 경계는 확장하지 못하여, 전체 분포역이 좁아지는 양상을 보인다(Kerr et al. 2015). 이는 온도에 민감한 한대성 곤충이 온난화에 취약함을 시사하며, 특히 분포의 남방 경계에 위치한 개체군은 서식에 적합한 기후조건을 상실할 위험이 크다(Hampe and Petit 2005).

우수리뒤영벌(Bombus ussurensis Radoszkowski, 1877)은 중국 동북부, 러시아 극동, 그리고 한반도에 분포하는 북방계 뒤영벌로, 남한은 이 종의 분포역 남방 경계에 해당한다. 국내 뒤영벌 분포 조사에 따르면, 우수리뒤영벌은 강원도 정선, 춘천, 평창 등 중부 및 북부 산간지역을 중심으로 채집되었으며 남부지역에서는 거의 발견되지 않아, 국내 여러 뒤영벌 중에서도 북방·산지에 편중된 분포를 보인다(Yoon et al. 2012). 또한 이 종은 전체 채집 개체 중 약 4.6%에 불과하여, 우점종인 B. ardens나 B. ignitus에 비해 개체군 규모가 작은 주변부 종으로 확인되었다. 따라서 우수리뒤영벌의 남한 개체군은 기후 온난화에 특히 취약할 것으로 예상된다. 최근 Rahimi and Jung (2024)이 국내 206종 화분매개곤충의 기후변화 분포를 예측하며 이 종을 포함한 바 있으나, 다수 분류군을 일괄 모델링한 광역 분석으로 우수리뒤영벌을 단일종 수준에서 심층적으로 다루거나 남단 경계 개체군의 취약성에 초점을 맞추지는 않았다. 더욱이 남한은 국토의 북쪽이 국경과 맞닿아 있어, 온난화에 따라 적합서식지가 북상하더라도 이를 수용할 수 있는 더 높은 위도의 지역이 국경 내에 존재하지 않는다는 지리적 특수성을 지닌다. 이러한 조건은 남단 개체군의 취약성을 더욱 심화시킬 수 있다.

종분포모형(Species Distribution Model, SDM)은 종의 출현자료와 환경변수의 관계를 통해 서식적합지를 추정하고 미래 기후조건에서의 분포 변화를 예측하는 도구로, 생태 및 보전 연구에 널리 활용되고 있다(Elith and Leathwick 2009). 그중 최대엔트로피 방법(MaxEnt)은 출현자료만을 이용하여 표본 수가 적은 경우에도 안정적인 예측 성능을 보여, 분포자료가 제한적인 종의 분석에 특히 유용하다(Phillips et al. 2006). 국내에서도 식물, 해양무척추동물 등 다양한 분류군을 대상으로 MaxEnt 기반 분포 예측 연구가 수행되어 왔으나(Park et al. 2018, 2020;Cho et al. 2020;Kim et al. 2022;Hong et al. 2023;Jeong et al. 2024), 화분매개곤충, 특히 기후변화에 취약한 한대성 뒤영벌을 대상으로 한 연구는 부족한 실정이다.

본 연구는 우수리뒤영벌을 대상으로 (1) 현재 기후조건에서의 남한 내 서식적합지를 규명하고, (2) 미래 기후변화 시나리오에 따른 서식적합지의 면적 변화와 분포 이동 방향을 예측하며, (3) 이를 바탕으로 남단 개체군의 보전학적 함의를 도출하고자 하였다. 이를 통해 기후변화에 따른 우수리뒤영벌 남한 개체군의 잠재적 분포 변화를 이해하고, 향후 분포 모니터링과 보전 전략 수립을 위한 기초자료를 제공하고자 한다.

2. 재료 및 방법

2.1. 종 출현 자료

본 연구의 대상종인 우수리뒤영벌(Bombus ussurensis Radoszkowski, 1877)은 중국 동북부, 러시아 극동, 한반도에 분포하는 북방계 뒤영벌로, 남한은 이 종 분포의 남방 경계에 해당한다. 종의 출현자료는 전지구생물다양성정보기구(Global Biodiversity Information Facility, GBIF)에서 남한지역의 지리좌표가 제공된 출현기록을 수집하였다(GBIF 2026). 수집된 원자료는 좌표 불확실성이 1 km를 초과하거나 육지 경계를 벗어난 기록을 제외하였으며, CoordinateCleaner 패키지(Zizka et al. 2019)의 cc_val(), cc_zero(), cc_dupl() 함수를 이용하여 좌표의 유효성을 검토하고 영점 좌표와 중복 좌표를 추가로 제거하였다. 정제된 자료 중 기후자료의 기준 기간(1981~2010)과 시간적으로 정합하도록 1980~2015년에 기록된 자료를 사용하였으며, 연도 정보가 없는 기록도 함께 포함하였다. CHELSA(Climatologies at High resolution for the Earth’s Land Surface Areas) 현재 기후자료가 30년 평균값이라는 점을 고려하여 기준 기간 전후 소수 연도의 기록을 포함하였고, 연도 미상 기록의 포함이 결과에 미치는 영향은 후술할 민감도 분석을 통해 검토하였다. 공간적 표본 편중을 완화하기 위해 spThin 패키지(Aiello-Lammens et al. 2015)를 이용하여 출현지점 간 최소거리 5 km로 설정한 공간 솎아내기(spatial thinning)를 수행하였다. 이를 100회 반복하여 가장 많은 출현지점이 유지된 자료세트를 선택하였으며, 최종적으로 92개 지점을 모델링에 사용하였다(Fig. 1). 최종 출현자료의 신뢰성을 확인하기 위해 자료의 출처를 검토한 결과, 최종 92개 지점 중 88개(95.7%)는 국립생태원(National Institute of Ecology)의 전문가 생태조사 자료였고, 나머지 4개(4.3%)는 보존표본 및 DNA 바코드 기반 기록이었다. 이들 4개 기록은 경기 및 충남지역에서 수집되었으며, 채집연도가 확인된 3개 기록은 각각 2005년, 2006년 및 2010년으로 현재 기후 기준기간(1981~2010)에 포함되었으며, 나머지 1개 기록은 채집연도가 확인되지 않았다. 아울러 표본·바코드 기반 기록과 관측기록이 여름철 최고기온 및 겨울철 최저기온으로 정의되는 기후공간에서 광범위하게 중첩됨을 확인하여 동정의 일관성을 검토하였다.

2.2. 환경변수

현재 및 미래 기후자료는 CHELSA version 2.1의 생물기후변수를 사용하였다(Karger et al. 2017). CHELSA는 약 1 km (30 arc-second) 공간해상도로 제공되며, 현재 기후는 1981~2010년 기준값이다. CHELSA는 지형물리 기반 다운스케일링을 적용하여 산악지형 등 지형이 복잡한 지역의 기후를 정밀하게 재현하고, 현재와 미래 자료가 동일한 다운스케일링 체계로 생성되어 시기 간 비교에 적합하다. 이러한 특성은 산간지역을 주 서식지로 하는 본 대상종의 분석에 적합하다. 초기 후보로 19개 생물기후변수 전체를 대상으로 출현지점에서의 값을 추출한 후, usdm 패키지(Naimi 2017)의 vifcor() 함수를 이용하여 변수 간 Pearson 상관계수(r)의 절대값 0.7을 기준으로 다중공선성을 검토하였다. 이후 대상종의 생태적 특성을 고려하여 최종 4개의 기후변수인 최난월 최고기온(bio5), 최한월 최저기온(bio6), 연강수량(bio12), 강수량 계절성(bio15)을 선정하였으며, 선정된 최종 변수들의 분산팽창계수(variance inflation factor, VIF)가 모두 10 미만임을 확인하였다.

서식지 조건을 반영하기 위해 ESA WorldCover 2021 version 2.0의 10 m 해상도 토지피복자료(Zanaga et al. 2022)로부터 세 가지 피복 비율 변수를 산출하였다. 각 토지피복 등급(산림, 초지, 시가지)에 대해 10 m 격자를 이진화한 후 1 km 격자 내 점유 비율로 집계하여 산림 비율(forest), 초지 비율(grassland), 시가지 비율(built-up)을 연속형 변수로 생성하고 기후변수와 동일한 격자에 정합하였다.

최종 예측변수는 기후 4종과 토지피복 3종을 합한 7개로 구성하였다(Table 1). 이때 토지피복변수는 시간에 따른 변화를 고려하지 않고 현재 상태로 고정하여 현재와 미래 예측에 동일하게 적용하였다.

2.3. 미래 기후 시나리오

미래 기후 조건은 CHELSA version 2.1이 제공하는 CMIP6 다운스케일 자료를 사용하였다(Karger et al. 2023). 이 자료는 CMIP6 기후모델 출력을 델타 변화법(delta-change method)에 기반하여 약 1 km 해상도로 다운스케일한 것이다. 기후모델 간 불확실성을 반영하기 위해 기후민감도가 상이한 3개 전지구기후모델(GFDL-ESM4, MPI-ESM1-2-HR, UKESM1-0-LL)을 선정하였다. 기후 시나리오로는 공통사회경제경로(Shared Socioeconomic Pathways, SSPs)인 SSP1-2.6과 SSP5-8.5를 적용하였다. SSP1-2.6은 지속가능한 발전과 친환경적 사회·경제 전환을 전제로 온실가스 배출이 낮게 유지되는 경로이며, SSP5-8.5는 화석연료에 의존하는 에너지 집약적 발전과 높은 온실가스 배출이 지속되는 경로이다. 두 시나리오는 2100년경 각각 약 2.6 W m^-2와 8.5 W m^-2의 복사강제력 수준에 해당한다. 미래 시기는 2050년대(2041~2070)와 2070년대(2071~2100)를 대상으로 하였다. 이에 따라 총 12개(3 GCM×2 SSP×2 시기)의 미래 기후 조합에 대해 예측을 수행하였다. 미래 예측에서는 토지이용이 현 수준으로 유지된다는 가정하에 현재의 토지피복변수를 적용하였다.

2.4. 종분포모형

종분포모형은 최대엔트로피 방법인 MaxEnt를 R 환경의 ENMeval 2.0 패키지(Kass et al. 2021)를 통해 maxnet 알고리즘으로 구축하였다. MaxEnt는 출현자료와 배경점의 환경조건 대비를 통해, 주어진 제약조건을 만족하는 분포 가운데 최대 엔트로피를 갖는 (즉 가장 균일한) 분포를 추정하는 기계학습 기법으로, 본 연구에서는 cloglog 출력을 서식적합도(0~1)로 사용하였다(Phillips et al. 2006;Elith and Leathwick 2009). 배경점(background points)은 대상종의 접근가능영역(accessible area, M)을 반영하기 위해 출현지점으로부터 100 km 완충영역 내에서 추출하였다(Barve et al. 2011). 이는 종이 여러 세대에 걸쳐 도달 가능한 범위로 배경점을 제한함으로써 모형이 서식 가능 공간을 적절히 대비하도록 한 것으로, 남한의 공간 규모를 고려한 설정이다. 이때 조사 편중을 보정하기 위해 GBIF의 벌목(Hymenoptera) 및 나비목(Lepidoptera) 출현기록의 공간밀도로 편중 레이어(bias layer)를 구축하고, 이에 비례하여 10,000개의 배경점을 추출하였다. 이는 특정 분류군의 생태적 유사성을 반영하기 위한 것이 아니라, 국내 곤충 조사 및 채집 노력의 공간적 편중을 대리하기 위한 표적집단 배경(target-group background) 접근에 따른 것이다(Phillips et al. 2009). 모형의 복잡도를 조절하고 과적합을 방지하기 위해 조절계수(regularization multiplier, RM)를 0.5~4.0 (0.5 간격)으로 설정하고, 특성계급(feature class)을 linear (L), quadratic (Q), hinge (H)를 조합한 L, LQ, LQH, H의 네 가지로 구성하였다. 변수 간 상호작용을 허용하여 모형의 복잡성을 증가시킬 수 있는 product (P)는 후보에서 제외하였으며, 불연속적인 반응형태를 생성하는 threshold (T) 역시 보다 단순하고 매끄러운 반응곡선을 유지하기 위해 제외하였다(Merow et al. 2013;Phillips et al. 2017). 조절계수와 특성계급의 조합으로 총 32개 후보 모형을 평가하였다. 공간적 자기상관을 고려하여 4분할 공간 블록 교차검증(spatial block cross-validation)을 적용하였으며, 최적 모형은 표본 수를 보정한 Akaike 정보기준(AICc)이 가장 낮은 모형을 우선 선정하고, 동일한 최소 AICc를 갖는 모형이 둘 이상인 경우 검증 AUC(Area Under the Curve)가 높은 모형을 선정하였다. 모형의 예측력은 검증 AUC, 누락률(omission rate), 연속 Boyce 지수(Continuous Boyce Index, CBI)로 평가하였다. 아울러 각 예측변수의 상대적 중요도는 최종 선정된 maxnet 모형에 대해 별도로 수행한 순열 기반 변수 중요도(permutation-based variable importance) 분석으로 평가하였다. 출현지점과 배경점의 환경자료를 결합한 평가자료에서 각 예측변수의 값을 하나씩 무작위로 치환한 후 cloglog 예측값을 다시 산출하고, 원래 예측값과 치환 후 예측값 간 Pearson 상관계수의 감소량(1-r)을 변수 중요도의 지표로 사용하였다. 각 변수에 대해 이 과정을 30회 반복하였으며, 음의 값은 0으로 처리한 후 평균하였다. 마지막으로 변수 간 상대적 중요도를 비교하기 위해 전체 변수의 중요도 합이 100%가 되도록 정규화하였다.

2.5. 서식적합지 분석

최종 모형으로부터 산출된 현재 및 미래의 서식적합도(cloglog-transformed suitability)에 대해 다음과 같은 분석을 수행하였다. 적합/비적합 구분에는 10 백분위수 훈련 존재(10th percentile training presence) 임계값을 적용하였으며, 분석에 사용한 임계값은 0.346이었다. 이를 기준으로 0.346 미만을 비적합으로 구분하고, 기존 종분포모형 연구의 적합도 분류방식(Thapa et al. 2018)을 참고하여 임계값 이상의 적합지역을 낮음(0.346 이상 0.50 미만), 중간(0.50 이상 0.75 미만), 높음(0.75 이상)의 세 등급으로 구분하였다. 등급별 면적 변화는 이들 4개 적합도 등급을 기준으로 산출하였으며, 면적가중 분포 중심과 평균고도 변화는 임계값(0.346) 이상의 적합지역을 대상으로 평가하였다.

미래 예측은 시나리오·시기별로 3개 GCM에 동일한 가중치를 부여하여 서식적합도 예측값을 픽셀 단위로 산술평균하고, 이를 앙상블 예측값으로 사용하였다. 동일 픽셀에서의 표준편차를 GCM 간 예측 불확실성의 지표로 산출하였다. 모형의 학습은 접근가능영역 전체에서 수행하되, 출현자료가 남한에 국한되므로 결과의 해석과 정량 분석은 남한지역으로 한정하였다. 이를 위해 GADM 행정경계 자료로 예측 결과를 마스킹하고, 위도별 격자 면적을 보정하여 실제 면적(km²)을 산출하였다. 현재 대비 미래의 적합지 변화는 유지, 축소, 확장으로 구분하였으며, 분포 이동의 방향성은 적합지의 면적가중 분포 중심(centroid)의 수평 이동거리·방위와 평균 고도 변화로 정량화하였다.

주 분석 결과가 분석 설정에 얼마나 민감한지를 평가하기 위한 보조 분석으로, 변수 선정 방식, 공간 솎아내기 거리 및 연도 미상 자료의 포함 여부를 달리하여 민감도 분석을 수행하였다. 각 대안 설정에서 모형을 재적합한 후 현재 서식적합도 예측값을 주 분석 결과와 픽셀 단위 Pearson 상관계수로 비교하였다. 민감도 분석의 세부 방법과 결과는 보충자료(Supplementary Materials)에 제시하였다. 모든 분석은 R version 4.5.1 (R Core Team 2025)에서 수행하였다.

3. 결 과

3.1. 모형 성능 및 변수 기여도

공간 블록 교차검증을 적용한 32개 후보 모형 중 AICc 기준으로 선정된 최적 모형은 힌지를 포함한 특성계급 조합(LQH)과 조절계수 2.0을 사용하였다(Table 2). 최적 모형의 검증 AUC는 0.805 (±0.110)로 양호한 수준이었으며, 10 백분위수 훈련 존재 누락률(OR10)은 0.022, 연속 Boyce 지수는 0.543이었다. 순열 중요도 분석 결과, 분포 결정에 가장 크게 기여한 변수는 최한월 최저기온(bio6, 45.1%)이었으며, 다음으로 산림 비율(forest, 29.7%), 강수량 계절성(bio15, 13.2%), 시가지 비율(built-up, 8.7%)의 순이었다(Table 3). 반응곡선상 최한월 최저기온은 약 -15°C 부근에서 서식적합도가 최대에 이른 후 온난한 조건에서 급격히 감소하였고, 산림 비율이 높을수록, 시가지 비율이 낮을수록 적합도가 증가하였다(Fig. 2). 최난월 최고기온(bio5)에 대한 반응은 약 25°C까지 증가한 후 안정되는 양상을 보였다.

3.2. 현재 서식적합지

현재 기후조건에서 비적합 지역은 52,708 km²로 남한 전체 분석면적의 49.9%를 차지하였으며, 적합지역은 낮음 20,230 km² (19.1%), 중간 20,109 km² (19.0%), 높음 12,672 km² (12.0%)로 나타났다(Table 4). 공간적으로 높은 적합도를 보이는 지역은 주로 강원도를 비롯한 동부 및 산악지역에 분포하였으며, 서부와 남부로 갈수록 적합도가 낮아지는 경향을 보였다(Fig. 3).

3.3. 미래 서식적합지 변화

미래 기후조건에서 전체 적합지역(낮음·중간·높음)의 면적은 모든 시나리오와 시기에서 현재보다 감소하였으며, 감소 폭은 후기 시기와 고배출 시나리오에서 더욱 크게 나타났다(Table 4, Fig. 4). 현재 대비 전체 적합지역의 감소율은 SSP1-2.6에서 2050년대에 36.7%, 2070년대에 45.5%였으며, SSP5-8.5에서는 각각 57.0%와 59.1%로 증가하였다. 적합도 등급별로도 높은 적합도를 보이는 지역의 감소가 뚜렷하였다. 고적합지역은 현재 12,672 km² (12.0%)에서 SSP1-2.6의 2050년대와 2070년대에 각각 8,554 km² (8.1%)와 7,824 km² (7.4%)로 감소하였으며, SSP5-8.5에서는 각각 6,745 km² (6.4%)와 5,207 km² (4.9%)로 감소하여, 고배출 시나리오의 후기에서 고적합지역의 축소가 가장 크게 나타났다.

3.4. 분포의 공간적 변화 및 이동

적합지역의 공간적 이동 역시 모든 미래 시나리오에서 일관된 방향성을 보였다(Table 5, Fig. 5). 면적가중 분포 중심은 현재 위치에서 북동쪽으로 이동하였으며, 이동거리는 SSP1-2.6의 2050년대와 2070년대에서 각각 18.7 km와 34.0 km, SSP5-8.5에서는 각각 41.6 km와 42.2 km로 나타났다. 적합지역의 평균고도 역시 모든 미래 조건에서 상승하였으며, 현재 대비 상승폭은 SSP1-2.6의 2050년대와 2070년대에서 각각 114 m와 147 m, SSP5-8.5에서는 각각 195 m와 206 m로 나타났다. 특히 고배출 시나리오에서 수평적 이동거리와 평균고도 상승폭이 모두 크게 나타났다.

3.5. 미래 분포 예측의 모델 간 불확실성

3개 GCM 간 서식적합도 예측의 표준편차는 지역과 시나리오에 따라 공간적으로 차이를 보였다(Fig. 6). 각 SSP에서 2070년대의 불확실성은 2050년대보다 전반적으로 넓은 지역에서 높게 나타났으며, 특히 SSP5-8.5의 2070년대에서 상대적으로 높은 GCM 간 변이가 가장 광범위하게 나타났다. 높은 불확실성은 주로 중·북부 및 산악지역을 중심으로 분포하였으며, 이는 미래 적합지의 세부적인 공간적 경계와 적합도 수준이 GCM에 따라 달라질 수 있음을 보여준다.

4. 고 찰

본 연구에서 예측된 Bombus ussurensis의 현재 서식적합지는 강원도 산간지역과 백두대간을 중심으로 형성되었으며, 이는 국내 뒤영벌 분포 조사에서 이 종이 정선, 춘천, 평창 등 중부 및 북부 산간을 중심으로 채집되고 남부지역에서는 거의 발견되지 않은 결과(Yoon et al. 2012)와 잘 일치하였다. 출현자료로 사용된 지점과 독립적인 현장 채집 기록이 모형 예측과 부합한다는 점은, 본 모형이 대상종의 실제 분포 경향을 신뢰성 있게 재현하였음을 뒷받침한다.

순열 중요도 분석에서 최한월 최저기온(bio6)이 가장 큰 기여를 보였고, 반응곡선상 약 -15°C 부근에서 서식적합도가 최대에 이른 후 온난한 조건에서 급격히 감소하였다. 이는 우수리뒤영벌이 충분히 한랭한 겨울 기후를 요구하는 북방계 종임을 정량적으로 보여준다. 뒤영벌은 한랭한 기후에 기원을 둔 분류군으로 온난화에 대한 취약성이 높으며, 온도가 이들의 지리적 분포를 결정하는 핵심 요인으로 알려져 있다(Kerr et al. 2015). 특히 뒤영벌은 여왕만이 월동(diapause)을 거쳐 이듬해 봄 단독으로 군체를 형성하므로(Yoon et al. 2012), 겨울철 기후조건은 개체군 유지에 직접적인 영향을 미치는 것으로 판단된다. 아울러 산림 비율이 높고 시가지 비율이 낮을수록 적합도가 증가한 결과는, 이 종이 산림 서식지를 선호하며 도시화된 경관을 회피함을 시사한다. 도시화와 집약적 토지이용에 따른 서식지 손실 및 밀원 감소는 뒤영벌 개체군을 위협하는 주요 요인으로 지목되어 왔다(Goulson et al. 2008;Cameron and Sadd 2020). 이러한 서식지 요건은 우수리뒤영벌이 산간지역에 편중되어 분포하는 경향(Yoon et al. 2012)과 부합한다.

미래 기후조건에서 우수리뒤영벌의 남한 내 서식적합지 면적은 현재(1981~2010) 대비 모든 시나리오에서 감소하였으며, 그 폭은 저배출 시나리오에서 36.7~45.5%, 고배출 시나리오에서 57.0~59.1%에 달하였다. 특히 높음 등급 적합지 비율이 현재 12.0%에서 SSP5-8.5 2070년대에 4.9%로 감소하고 비적합지가 49.9%에서 79.5%로 증가하는 등, 전 등급에 걸쳐 단계적 서식지 질 저하가 관찰되었다. 배출 시나리오가 높고 시기가 늦을수록 감소 폭이 커지는 일관된 경향은 우수리뒤영벌의 분포가 온난화 강도에 직접적으로 반응함을 보여준다. 이는 기후 온난화에 따른 뒤영벌의 분포역이 남쪽 경계에서 좁아진다는 대륙 규모의 관찰(Kerr et al. 2015)과 궤를 같이하며, 국내 남단 개체군에서도 동일한 취약성이 나타남을 확인한 결과이다. 본 연구의 예측은 국내 화분매개곤충을 대상으로 한 Rahimi and Jung (2024)의 결과와 상이하였다. Rahimi and Jung (2024)은 우수리뒤영벌의 미래 적합지역이 기후변화에 따라 확대되고 남쪽으로 이동할 것으로 예측한 반면, 본 연구에서는 모든 시나리오에서 적합지역이 감소하고 잔존 적합지가 북동쪽 및 고지대로 이동하는 것으로 나타났다. 이러한 상반된 예측에는 두 연구 간 예측변수 구성과 기후자료, 출현자료 및 모형 설정의 차이가 복합적으로 영향을 미쳤을 가능성이 있다. 특히 Rahimi and Jung (2024)은 등온성, 연교차, 최난분기 평균기온 및 강수 관련 변수를 사용하였으나 최한월 최저기온(bio6)은 포함하지 않았다. 반면 본 연구에서는 bio6이 높은 변수 중요도를 보였으며, 반응곡선에서도 적합도와 뚜렷한 비선형 관계를 나타내어 우수리뒤영벌의 잠재분포를 설명하는 주요 기후변수로 확인되었다. 따라서 겨울철 기온을 포함한 예측변수 구성의 차이가 두 연구의 상반된 미래 예측에 기여했을 가능성이 있다. 또한 두 연구는 사용한 기후자료(WorldClim 대 CHELSA), 출현자료의 구성 및 분석 범위(다분류군 분석 대 단일종 분석)에서도 차이가 있으므로, 예측 결과의 차이를 특정 변수 하나의 영향으로 단정하기는 어렵다. 이러한 결과는 기후변화에 따른 종분포 예측이 환경변수의 선정과 모형 설정에 민감할 수 있으며, 대상종의 생태적 특성을 고려한 변수 선정과 예측 불확실성의 검토가 중요함을 보여준다. 실제로 본 연구의 민감도 분석에서도 공간 솎아내기 거리와 연도 미상 자료의 처리에는 예측이 비교적 강건한 반면, 예측변수의 선정 방식에 따라 공간적 예측이 크게 달라져(Supplementary Table S1, Fig. S1), 대상종의 생태적 특성을 반영한 변수 선정이 중요함을 뒷받침하였다.

본 연구의 가장 주목할 만한 결과는, 광역 예측 범위에서는 적합지가 북한 북부 및 그 이북의 고위도지역으로 이동하였으나 남한 국경 내에서는 이러한 북쪽으로 이동을 수용한 지역이 사실상 부재하였다는 점이다(Supplementary Figs. S2 and S3). 남한 내 신규 확장 면적은 모든 시나리오에서 매우 제한적(157~214 km²)이었으며, 이에 따라 축소가 확장을 압도하여 서식지의 순손실로 귀결되었다. 이는 온난화에 따라 적합서식지가 북상하면서 그 상당 부분이 남한의 북쪽 국경을 넘어 형성되기 때문으로, 그 결과 남한 영역 내에서는 새로이 형성되는 적합지가 거의 없이 기존 서식지의 손실만 나타나게 된다. 즉, 종의 잠재적 분포 자체는 북쪽으로 이동하나, 남한이라는 행정적 단위를 기준으로 보면 이러한 북상이 국경 외부에서 일어나 국내 서식지의 순손실로 귀결되는 것으로 해석된다.

한편 적합지의 면적가중 분포 중심은 미래에 북동 방향으로 18.7~42.2 km 이동하였고, 평균 고도는 현재 404 m에서 114~206 m 상승하였다. 이는 수평적 북상 경로가 차단된 상황에서 잔존 적합지가 보다 서늘한 고지대로 집중되는 수직적 이동(upslope shift)이 진행됨을 의미한다. 뒤영벌이 온난화에 대응하여 고도가 높은 지역으로 이동한다는 것은 대륙 규모 분석에서도 확인된 바 있다(Kerr et al. 2015). 그러나 산악 지형에서는 고도가 높아질수록 이용 가능한 면적이 급격히 감소하므로(Elsen and Tingley 2015), 고도 상승에 의한 피난은 결국 서식지가 산 정상부에 고립되어 소멸에 이르는 이른바 ‘멸종으로의 에스컬레이터(escalator to extinction)’ 현상으로 이어질 수 있다(Freeman et al. 2018). 수평·수직 양방향에서 피난처가 제약되는 남한 개체군의 상황은, 이 종의 국지적 취약성이 여타 지역보다 클 수 있음을 시사한다.

본 연구의 예측은 여러 측면에서 보수적으로 해석될 필요가 있다. 또한 미래 서식적합도 예측에는 GCM 간 불확실성이 존재하였으며, 그 정도는 시나리오와 시기에 따라 공간적으로 달랐다. 특히 후기 시기에는 일부 중·북부 및 산악지역에서 상대적으로 높은 GCM 간 변이가 나타났으므로(Fig. 6), 미래 적합지의 세부적인 공간 경계는 단일한 확정적 예측이라기보다 GCM 간 불확실성을 고려하여 해석할 필요가 있다. 반응곡선상 최난월 최고기온(bio5)에 대한 적합도는 약 25°C 이후 감소하지 않고 안정되는 양상을 보였는데, 이는 남한이 이 종 분포의 남방 경계에 해당하여 출현자료가 종의 고온 내성 상한을 충분히 포괄하지 못하는 니치 절단(niche truncation)에 기인한 것으로 판단된다(Chevalier et al. 2022). 지역 규모 자료로 구축된 종분포모형은 현재 조건에 대한 예측에서는 신뢰할 만하나, 미래 기후와 같이 학습 자료의 범위를 벗어나는 조건으로 외삽할 경우 서식지 변화를 실제와 다르게 추정할 수 있다(Chevalier et al. 2022). 또한 본 연구에서는 미래 기후조건의 변화에 따른 잠재분포 변화를 평가하는 데 초점을 두어 토지피복은 현재 상태로 고정하였다. 그러나 SSP 기반의 미래 토지이용자료가 이용 가능하며, 토지이용 변화는 기후변화와 상호작용하여 뒤영벌의 미래 분포에 영향을 미칠 수 있다(Marshall et al. 2018). 따라서 미래 토지피복의 변화를 반영하지 못한 것은 본 연구의 중요한 한계이며, 본 연구의 미래 예측은 토지피복이 현재 상태로 유지된다는 가정하에서 나타나는 기후변화에 따른 잠재적 분포 변화로 해석할 필요가 있다. 향후에는 SSP 기반 미래 토지이용 시나리오를 기후 시나리오와 함께 적용하여 기후와 토지이용 변화의 개별 및 복합적인 영향을 평가할 필요가 있다. 또한, 조사 편중 보정을 위해 사용한 GBIF 벌목, 나비목 자료는 조사시기, 채집방법 및 조사기관이 서로 다를 수 있어 실제 조사노력의 공간적 편중을 완전히 반영하지 못했을 가능성이 있으며, 향후에는 대상종과 조사시기 및 채집특성이 유사한 분류군의 자료를 활용한 검증이 필요하다.

본 연구는 기후와 토지피복을 중심으로 서식적합성을 평가하였으므로, 뒤영벌의 분포를 실제로 제약하는 밀원식물의 가용성, 둥지 자원, 종간 상호작용 등 생물적 요인을 직접 반영하지 못하였다. 아울러 출현자료의 수집 기간(1980~2015년)이 기후 기준 기간(1981~2010년)과 완전히 일치하지는 않으며 일부 연도 미상 자료를 포함하였으나, 연도 미상 자료의 포함 여부에 따른 민감도 분석에서도 현재 서식적합 예측은 주 분석 결과와 공간적 상관(r=0.96)을 보여(Supplementary Table S1, Fig. S1), 이로 인한 영향은 제한적인 것으로 판단된다. 나아가 출현자료가 남한에 국한되어 종 전체의 기후 니치를 완전히 포괄하지 못한 한계가 있으며, 향후 중국 동북부 및 러시아 극동을 포함한 전 분포역 자료를 활용한 분석이 이루어진다면 예측의 정확성을 높일 수 있을 것이다.

본 연구 결과는 우수리뒤영벌의 남한 개체군이 기후변화에 매우 취약하며 그 위협이 금세기 내에 현실화될 수 있음을 보여준다. 이 종의 잔존 적합지가 강원 산간 고지대에 집중될 것으로 예측되므로, 향후 보전 논의에서 해당 지역이 우선적으로 고려될 필요가 있다. 다만 현재 우수리뒤영벌은 국내에서 보호종으로 지정되어 있지 않고 이 종에 특화된 분포 및 개체군 자료도 부족한 실정이므로, 우선적으로는 강원 산간지역을 중심으로 한 분포 실태조사와 장기 모니터링을 통해 기초자료를 축적하는 것이 필요하다. 본 연구에서 제시한 우수리뒤영벌의 현재 및 미래 서식적합지 예측은 기후변화에 취약한 지역을 파악하고, 향후 모니터링 우선 대상지 설정과 보전 전략 수립을 위한 기초자료로 활용될 수 있다.

적 요

기후변화는 한대성 화분매개곤충의 서식지에 큰 영향을 미치며, 분포역 남방 경계에 위치한 개체군은 특히 취약하다. 본 연구는 중국 동북부와 러시아 극동, 한반도에 분포하는 북방계 뒤영벌 우수리뒤영벌(Bombus ussurensis)을 대상으로, 현재 및 미래 기후조건에서의 남한 내 서식적합지 변화와 분포 이동을 예측하였다. 최대엔트로피 모형(MaxEnt)을 이용하여 기후변수 4종과 토지피복변수 3종을 바탕으로 서식적합지를 추정하였으며, 3개 전 지구 기후모델과 2개 공통사회경제경로(SSP1-2.6, SSP5-8.5), 2개 미래 시기(2050년대, 2070년대)에 대해 투영하였다. 모형에서 최한월 최저기온이 가장 중요한 환경요인으로 나타나, 이 종이 한랭한 겨울 기후를 요구하는 북방계 종임을 확인하였다. 현재 서식적합지는 강원 산간지역을 중심으로 분포하였으며, 이는 기존 현장 조사 결과와 일치하였다. 미래에는 모든 시나리오에서 서식적합지가 감소하였고(현재 대비 36.7~59.1%), 고배출 시나리오와 장기 시기일수록 감소 폭이 컸다. 적합지의 면적가중 분포 중심은 북동 방향과 고지대로 이동하였으나, 적합서식지의 북상이 상당 부분 남한의 북쪽 국경 밖에서 일어나 남한 영역 내에서는 서식지 확장이 거의 이루어지지 못하였다. 이러한 결과는 남한의 우수리뒤영벌 개체군이 수평·수직 양방향에서 기후 피난처가 제약된 상황에 놓여 있어 기후변화에 매우 취약함을 시사한다. 본 연구는 우수리뒤영벌의 현재 및 미래 잠재분포 변화를 정량적으로 제시함으로써, 기후변화에 따른 화분매개곤충의 취약성을 이해하고 향후 보전 전략을 수립하기 위한 기초자료를 제공한다.