서 론
너도밤나무속(Fagus ssp.)은 참나무과(fagaceae) 낙엽활 엽교목으로 전세계 약 10~11종이 분포하고 있으며(Denk et al., 2005;Fang and Lechowicz, 2006), 북반구 온대 활엽 수림의 대표적인 속(genus) 중 하나이다(Shen, 1992;Denk, 2003). 우리나라와 인접한 너도밤나무속으로 일본산의 Fagus crenata Blume, Fagus japonica Maxim., 중국산의 Fagus engleriana Seemen ex Diels가 분포하고 있다. 우리나라의 경우 울릉도에만 한정적으로 너도밤나무(Fagus multinervis Nakai)가 자생하고 있으며, 산림청 희귀식물목록(2012년) LC등급(약관심종, Least Converned)으로 분류되고 있다.
울릉도 너도밤나무 생장 특성으로 Lee et al.(2007)은 울릉도에 우점하고 있는 대부분의 너도밤나무는 맹아로 갱신되며 큰 개체로부터 공급된 종자 발아는 광요구도가 낮아 성숙목이 있어도 생존이 가능하나, 성장함에 따라 광요구도가 높아지며 외적인 방해가 없는 한 요구조건을 충족시키기 어려워 생장을 멈추거나 고사한다고 알려져 있다. 너도밤나무속과 관련한 기후적인 특징으로 Yim(1983)에 따르면 북반구의 너도밤나무는 연평균 기온 4.5°C ~ 20.0°C, 뚜렷한 건기가 없는 선에서 연평균 강수량 600 ~ 1,000㎜ 범위에서 발견된다고 보고 있다. Fang and Lechowicz(2006)에서는 너도밤나무속의 지리적 분포 한계를 결정짓는 요인이 온도와 수분이며, 일부 종의 경우에만 역사적 요인이 결정한다고 주장한다. 너도밤나무가 자생하는 울릉도의 기후는 대륙성 기후인 한반도 내륙과 달리 해양성 기후로서 연중 온난하고 강수량 또한 사계절 고루 분포되어 동일 위도상의 다른 지역에 비하여 온난한 특징을 가지고 있다(Oh, 1978;Lee et al., 2007).
최근 포항시 금광리 일대에 여러 식물 화석 중 너도밤나무가 발견되었는데, 과거 한반도 육지에도 너도밤나무가 생육하였음을 추측할 수 있다(Paik et al., 2010). 한국의 너도밤나무와 관련한 선행연구로는 울릉도 너도밤나무림의 생태적 특성에 관한 연구(Han et al., 2019), 울릉도 군집구조와 하층식생 분포특성(Cheon et al., 2012), 울릉도 너도밤나무 하위군락별 토양특성(Park et al., 2000) 등이 존재한다. 울릉도를 대상지로 한 연구는 진행되었으나 울릉도에만 분포하고 있는 너도밤나무에 대하여 육지에서의 생육 가능성과 관련한 연구는 부재하다.
반면 우리나라와 인접한 일본은 혼슈 북부, 시코쿠, 규슈의 태평양 쪽에 F. japonica가 분포하고 있으며, 홋카이도 남부부터 규슈까지 광범위하게 F. crenata가 분포하고 있다 (Hara, 2010;Fang and Lechowicz, 2006). 관련한 선행연구로 기후변화에 따른 F. crenata의 민감성과 취약성 연구 (Matsui et al., 2004a), 기후 통제에 따른 F. crenata의 분포 연구(Matsui et al., 2004b) 등 종에 대한 분포, 특성과 관련한 연구가 진행되었다. 이외에도 F. grandifolia 분포 예측을 위한 종분포모형(Flessner et al., 2017), 기후변화에 따른 F. orientalis의 과거, 현재, 미래 분포 예측 연구(Dagtekin et al., 2020), 화석 기록과 전통적인 판별 지표를 활용한 F. sylvatica 종분포모형(Poli et al., 2022) 등 여러 나라에서 너도밤나무속을 대상으로 분포 확률, 기후변화 시나리오에 따른 분포 변화 등 연구된 바 있다.
종 출현 가능성을 도출하는 모형인 종분포모형(Species Distribution Model, SDM)은 생물종의 출현·비출현 데이터와 환경 데이터 사이의 관계를 분석한다. 생물 종 분포에 대한 정량적인 예측이 가능하여 최근에는 인간이 접근할 수 없는 접경지역의 생물종 예측, 재배식물·외래식물 예측, 기후변화 시나리오에 따른 생물종 예측, 산사태 발생 가능성 등 다양한 연구에 활용되고 있다(Sung et al., 2018;Lee et al., 2012;Lee et al., 2016;Kim et al., 2013).
한편 식물의 지리적 분포에는 기후적 인자와 지형, 토양과 같은 비기후적 인자 등 다양한 요소가 관여하는데(Chung, 2013), 가장 중요한 요인은 기후적 요인이다(Woodward, 1987;Woodward and Mckee, 1991). 이에 본 연구에서는 식물의 지리적 분포 결정에 우선이 되는 기후적인 관점으로 우리나라 육지에서 너도밤나무속의 생육 가능 확률과 생육 환경 특성을 살펴보고자 한다.
연구방법
1. 연구범위 및 연구대상지
너도밤나무속(Fagus spp.)에 대한 분석환경과 종의 생태적 특성을 고려하여 연구범위와 연구대상지를 설정하였다. Oh(2015)의 연구결과에 따르면 F. multinervis와 계통이 가까운 F. engleriana, F. crenata, F. japonica 3종과 유전자 계통분석을 수행한 결과, F. multinervis와 F. japonica, F. engleriana 가 자매관계를 이루는 것으로 나타났다. 이를 통해 서로가 생태적 특성을 공유할 수 있는 가능성을 두고 연구대상종를 설정하였다. 다만, F. engleriana의 경우, 출현 데이터 확보가 어려움에 따라 Table 1과 같이 유형을 분류하였다. 연구대상 종은 울릉도에 생육하는 F. multinervis, 일본에 생육하는 F. crenata, F. japonica로 3종을 선정하였다. 우리나라에 자생하는 F. multinervis만을 분석한 경우(Case1_Fm), F. multinervis와 자매관계를 이루는 F. japonica를 한 종으로 묶어 분석한 경우(Case2_Fmj), F. multinervis와 근연종이 아닌 F. crenata를 묶어 분석한 경우(Case3_Fmc), F. multinervis, F. japonica, F. crenata 세 종을 묶어 너도밤나무속으로 분석한 경우(Case4_Fmjc), 3종 각각을 분석한 뒤 중첩한 경우(Case5_Ovlp) 까지 총 5종류이다. 연구대상지는 유형에 따라 Case1_Fm은 한반도 남한, Case2_Fmj~Case5_Ovlp는 한반도 남한, 일본으로 설정하였다.
2. 종분포모형 분석
1) 종분포모형
본 연구에서는 우리나라 육지에서 너도밤나무 생육이 가능한 지역을 알아보기 위해 종분포모형의 다양한 모형 중 Maxent(maximum entropy) 모형을 이용하였다. Maxent 모형은 최대 엔트로피에 기반한 종분포모형으로 일반적인 종의 출현·비출현 데이터를 이용해야 하는 다른 모형과 달리 출현 데이터만으로도 예측의 정확도가 높은 기계학습 (machine learning) 방법 중 하나이다(Elith et al., 2006;Gibson et al., 2007;Seo et al., 2008).
예측 결과의 정확도를 높이기 위해 학습 데이터와 검증 데이터를 분할하는 교차검증(cross-validation)을 10번 반복 하였다. 교차검증은 모형의 정확도는 수신자 조작 특성(Recevier Operating Characteristic, ROC) 곡선의 하부면적 값(Area Under Cover, AUC)을 이용해 측정하였다. AUC값을 이용한 정확도는 기준값에 독립적인 장점을 가지고 있으며, 다양한 모형의 정확도에 많이 이용된다(Thuiller, 2003;Franklin, 2009). Maxent 알고리즘에 의한 기여도 평가를 살펴보기 위해 백분율 기여도(percent contribution)와 무작위 중요도(permutation importance)를 살펴보았다. 기여도 평가는 모형이 실측자료를 이용한 학습 중 각 변수가 결과 도출에 기여하는 정도를 의미하며, 모형 알고리즘에 따라 달라질 수 있는 상대적인 기여도이다(Lee et al., 2015). 이후 생육 확률에 대한 패턴을 파악하기 위해 확률을 서식지 부적합 지역(0-20%), 적합성이 낮은 지역(20-40%), 적합성이 중간인 지역(40-60%), 적합성이 높은 지역 (60-80%), 아주 적합한 지역(80-100%)으로 5개 범주로 분류하였다(Sim et al., 2022;Lee et al., 2021, Table 2). 분석에 대한 전체 흐름은 Figure 1과 같다.
2) 종 출현 데이터
너도밤나무속 출현자료는 세계생물다양성정보기구(Global Biodiversity Information Facility, GBIF)를 통해 F. multinervis, F. japonica, F. crenata 3종에 대한 위·경도 데이터를 확보하였다. 인간 관측(human observation)만으로 GBIF 데이터를 확보한 후 불명확한 출현 지점, 결측값, 중복된 좌표를 제거하였다. 이후 표본 편향(samplingbias)에 의한 공간적 자기상관(spatial autocorrelation)이 발생할 수 있으므로 한 픽셀에 하나의 위치 데이터가 들어가도록 중복된 지점을 제거하여 최종 출현 자료를 선정하였다 (Betts et al., 2006;Segurado et al., 2006). 최종 분석 데이터로 F. multinervis 14개, F. japonica 206개, F. crenata 446개를 선정하였다. 출현 데이터 전처리 과정은 R 3.3.2 프로그램을 사용하였다.
3) 환경 데이터
생육 확률을 분석하기 위한 환경변수로 Worldclim(https:// www.worldclim.org/)에서 제공하는 1970-2000년까지의 생물기후데이터 Bioclim (bioclimatic variables) 자료를 내려받아 활용하였다. Bioclim은 온도 및 강수량 데이터로부터 도출된 19개의 기후자료이며, 30초(약 1㎢)로 공간해상도를 가진다(Fick and Hijmans, 2017;Harris et al., 2020). 19개 변수 중 상관관계가 높은 변수들을 포함하면 다중공선성(collinearity)이 발생할 수 있으며 결과가 왜곡될 수 있다 (Watling et al., 2012). 상관성이 높은 값들을 제거하기 위해 분산팽창계수(Variance Inflation Factor, VIF)를 이용하였다(Kwon and Kim, 2023). 분산팽창계수는 다중공선성을 탐색하고 변수를 선택하기 위한 통계적인 지표로 일반적으로 VIF 10을 기준으로 판단한다(Park and Sung, 2010). 상관성이 높은 값들을 제거하여 최종변수 8개(Bio3, Bio5, Bio7, Bio8, Bio9. Bio14, Bio15, Bio18)를 선정하였다 (Table 3). 환경 데이터 편집, VIF 분석은 QGIS 3.34.8, R 3.3.2 프로그램을 사용하였다.
3. 생육환경 특성 분석
기후특성을 살펴보기 위해 생물기후데이터 Bioclim을 사 용하였다. 19개 변수 중 종분포모형 분석에 사용된 환경변수 Bio3, Bio5, Bio7, Bio8, Bio9, Bio14, Bio15, Bio18과 기후의 기초적인 지표인 연평균기온(Bio1), 연간 강수량(Bio12)까지 총 변수 10개를 선택하였다. 분석 대상은 6곳으로 우리나라 전체, 우리나라 육지에서 생육 가능성이 나타난 곳에 대하여 bioclim 래스터를 추출하여 각 변수마다의 평균을 구하였다. 또한, GBIF에서 획득한 F. multinervis, F. japonica, F. crenata 각 종마다의 위·경도 출현지점, 3종을 모두 합친 출현지점까지 총 4곳의 Bioclim 값을 추출하여(point sampling) 평균을 계산하였다. 픽셀 값 추출은 QGIS 3.34.8 Point sampling tool을 활용하였다.
결과 및 고찰
1. 생육가능지역 분석
1) 정확도 및 기여도 평가
종분포모형의 정확도 ROC의 AUC값은 Case1_Fm 0.946, Case2_Fmj 0.923, Case3_Fmc 0.876, Case4_Fmjc 0.871로 유형 모두가 신뢰할 수 있는 수준으로 나타났다(Figure 2). AUC값은 최소 0.5를 기준으로 하며, 1에 가까울수록 모형의 예측 정확도가 높으며, 일반적으로 0.8 이상이면 모형 예측력이 우수하다고 평가한다(Thuiller, 2003;Franklin, 2009).
해당 모형에 어떠한 변수가 영향을 미쳤는지에 대한 변수기여도를 나타내었다(Table 4). 백분율 기여도(percent contribution)는 모델 학습과정에서 변수 기여도로 각 환경 변수가 Maxent 모형의 적합도에 얼마나 기여하는지를 나타낸다. 무작위 중요도(permutation importance)는 입력변수 중 하나의 영향을 제거하여 성능의 차이가 얼마나 나는지 확인하는 방법(Lee et al., 2012)으로 훈련 데이터에서 변수의 값들을 섞어서 AUC의 감소를 측정하고, 큰 감소가 발생하면 해당 변수가 모델에 중요하다는 것을 나타낸다 (Phillips et al., 2006). Case1_Fm의 경우 종분포모형 분석 결과 우리나라 육지에서 생육 가능성이 거의 관측되지 않았으며, Case5_Ovlp의 경우 각 종마다의 분석 결과를 중첩한 것이기 때문에 변수기여도를 나타내지 않았다. Case1_Fm, Case5_Ovlp을 제외한 Case2_Fmj, Case3_Fmc, Case4_Fmjc 에서의 백분율 기여도를 살펴보면, Case2_Fmj의 경우 Bio3, Bio7, Bio5 순으로, Case3_Fmc의 경우 Bio5, Bio18, Bio15 순으로, Case4_Fmjc의 경우 Bio5, Bio18, Bio15 순으로 영향을 높게 받았다. 무작위 중요도는 Case2_Fmj에서 Bio15, Bio5, Bio7 순으로, Case3_Fmc에서 Bio5, Bio15, Bio18 순으로, Case4_Fmjc에서 Bio15, Bio5, Bio18 순으로 기여도가 높았다. 너도밤나무속의 생육 확률 분포를 결정하는 환경변수는 공통적으로 가장 따뜻한 분기의 최고기온(Bio5), 강수량 계절성(Bio15)의 기여도가 높았다.
2) 생육 가능성 예측 결과
생육 가능성 예측 결과는 Figure 3과 같다. 모형의 예측 결과는 0에서 1 사이의 확률 분포 지도로 나타나며, 1에 가까울수록 생육 확률이 높음을 의미한다. 우리나라만을 대상으로 한 Case1_Fm에서는 우리나라에서 울릉도를 제외하고 생육 가능성이 거의 나타나지 않았다. 이는 너도밤나무 출현 데이터가 공간적 자기상관 제거 이후 종분포모형 분석을 위한 최소 개체수만큼 적으며 비교적 좁은 지역에 집중되어 있기 때문이라 추측된다. 멸종위기종 종분포모형의 경우 분석을 위해서는 최소 출현지점이 7개 이상이 되어야 유의미한 결과가 나오는 것과 유사한 사례라 판단된다 (Franklin, 2009;Kim et al., 2012). Case2_Fmj, Case3_Fmc 에서 F. japonica, F. crenata에 대한 일본 생육 예측 결과 정확도를 살펴보기 위해, 식물사회학적 식생자료인 PRDB (Phytosociological Relevé Database)에서 제공하는 식물분 포도(PRDB map)와 비교하였다. 비교 결과, Case2_Fmj 생육 확률과 F. japonica의 분포, Case3_Fmc 생육 확률과 F. crenata 분포가 비슷한 결과로 나타났다(Tanaka, 2007). Case2_Fmj는 삼척시, 울진 등 동해안 가장자리에 생육 가능성이 발견되었으며, Case3_Fmc은 타 유형에 비해 동해안 지역에서 생육 가능성이 낮게 관측되었다(Figure 4). Case4_Fmjc는 Case2_Fmj의 F. japonica, Case3_Fmc의 F. crenata 두 종에 대한 예측 확률이 겹쳐 도출되었다. Case4_Fmjc에서 생육 가능성은 Case2_Fmj와 Case3_Fmc 의 중간 정도로 동해안 가장자리에 낮은 생육 가능성이 관측되었다. 마지막으로 Case5_Ovlp은 타 유형과 동일한 경향으로 동해안 지역에서 낮은 생육 가능성이 관측되었다.
2. 너도밤나무 분포 환경
너도밤나무속에 대한 기후 특성을 파악하기 위해 6가지로 구분하였다. 우리나라 전체 지역, 우리나라 육지에서 생육 가능성이 나타난 동해안 지역, 마지막으로 F. multinervis, F. japonica, F. crenata 3종에 대한 각각의 실제 출현지점, 3종 전체까지 총 6개 분류에 대해 평균과 표준편차를 계산 하였다(Table 5). 변수기여도가 높았던 Bio15(강수량 계절성)의 경우 강수량이 얼마나 변동되는가에 대한 지표로서 값이 클수록 변동이 높음을 의미한다. 우리나라, 동해안 지역, F. multinervis를 살펴보면, F. multinervis가 변동이 가장 낮았으며, 동해안 지역, 우리나라 순으로 변동이 높아졌다. 기여도가 높은 또다른 환경변수 Bio5(가장 따뜻한 달의 최고기온)는 우리나라가 가장 높았으며 다음으로 동해안 지역, F. multinervis순으로 낮았다. F. multinervis에 대한 기후 특성을 살펴본 결과, Bio1(연평균 기온)은 11.06±0.91, Bio12(연간 강수량)은 1299.95±48.11로 나타 났다. F. multinervis, F. japonica, F. crenata 세 종에 대한 환경 특성을 살펴본 결과, F. multinervis가 연간 강수량 (Bio12), 강수량 계절성(Bio15)의 값이 가장 낮았다. 동해안 지역의 기후를 살펴보면 Bio7(연간 기온 범위)가 우리나라 전체 평균보다 낮았는데, 기온 변동성이 우리나라의 다른 곳과 달리 동해안지역이 계절적인 온도 차가 평균적으로 낮음을 보여준다. 또한, 우리나라 평균에 비해 가장 건조한 분기의 평균기온(Bio9)과 최고 건조한 달의 강수량(Bio14)은 높고, 강수량 계절성(Bio15)은 낮았다. 이를 통해, 우리나라 육지에서 생육 가능성이 관측된 동해안 지역은 우리나라 전체에 비해 계절적인 온도 차가 낮고, 강수량 변동이 높지 않아 계절성이 강하지 않은 걸로 판단된다.
3. 종합고찰
본 연구는 우리나라 육지에서 너도밤나무속의 생육 가능성을 파악하기 위해 Maxent 모형을 활용하여 생육 확률을 분석하였다. 우리나라의 경우 울릉도에만 너도밤나무가 분포하기 때문에 우리나라 인근의 너도밤나무속 중 자료 수집 이 가능했던 F. japonica, F. crenata 출현 종과 함께 5가지로 분류하여 분석하였다. 분석 결과, Case1_Fm은 육지에서의 생육 가능성이 매우 낮았다. Case2_Fmj~Case5_Ovlp은 공통적으로 육지에서의 생육 가능성은 낮았으나 동해안 지역에서 동일하게 생육 가능성이 일부 관찰되는 경향이 보였다. 특히, 우리나라 너도밤나무와 자매관계인 F. japonica와 함께 분석하였을 때(Case2_Fmj), 타 유형에 비해 동해안에 서의 생육 가능성이 비교적 높게 나타났다. 이를 기후적인 측면에서 종합해보면 우리나라 육지에서는 전반적으로 너도밤나무 생육이 어렵지만 동해안 일부 지역에서 생육이 가능한 것으로 판단된다.
Case2_Fmj~Case5_Ovlp에 대한 울릉도 생육 확률을 살펴보았을 때 출현 지점이 존재함에도 불구하고 평균 생육 확률이 서식지 적합성이 높은 지역, 아주 적합한 지역이 적게 나타났다. 이는 환경변수의 패턴 학습, 샘플링 편향, 모델의 과적합 등 다양한 근거가 존재할 수 있다. 이번 분석의 경우 울릉도의 출현 지점에의 환경변수 래스터에 비해 일본의 환경변수 래스터가 주축으로 되어 모델이 학습된 것이라 판단된다.
너도밤나무속 3종에 대한 분포 환경 특성을 분석하였다. 다만, 집단마다 표본 개체수가 상이하므로 정밀한 통계 분석이 필요하다. 또한, 본 연구에서는 Worldclim의 Bioclim 데이터를 활용하여 기후관점으로 분포 상황을 살펴보았으나 추후 우리나라 기상청(KMA; Korea Meteorological Administration) 방재기상관측장비(AWS; Automatic Weather System)와 종관기상관측장비(ASOS; Automated Synoptic Observing System)에서 측정한 기후자료, 일본 기상청 (Japan Meteorological Agency, JMA) 자료를 활용하여 실제 기상 데이터 비교를 통한 정밀한 분석이 필요하다.
이번 연구는 기후적인 관점으로 보았기 때문에 Bioclim 데이터를 사용하였으나, 추후 우리나라 생육 특성을 반영할 수 있는 고해상도 기후변수와 함께 지형, 토양 등의 다양한 변수가 투입되어야 할 것이다. 종과 관련하여도 우리나라 너도밤나무 1종, 일본산 2종으로 총 3종에 대해 진행하였으 나 F. engleriana에 대한 추가적인 데이터를 획득하여 동아 시아 너도밤나무속에 대한 정리가 필요해 보인다. 또한 종 분포모형 중 종의 출현 데이터만을 이용하여 분석하는 Maxent 모형을 사용하였다. 최근에는 Maxent 모형 이외에도 여러 분석기법이 사용되고 있는데, 특히 여러 모형들을 분석하여 최적의 분석 결과를 도출하는 앙상블(ensemble) 기법이 많이 활용되고 있는 추세이다(Kwon, 2014;Chung et al., 2020). 추후 앙상블 기법을 활용하여 분석이 진행된다면 너도밤나무속에 대한 더욱 정확한 결과가 도출될 것이라 기대한다.
이렇듯 너도밤나무속에 대한 새로운 생육가능지역 발굴에 대한 연구는 또다른 분포 모형 연구로 연결될 것이다. 생물종 분포의 예측 연구는 대상이 자연 현상이기 때문에 모형의 수립과 적용의 결과가 후속 모형 연구의 방향을 위한 토대가 된다(Kim et al., 2015). 때문에 종분포모형에 사용되는 환경요인, 분석기법 등에 대해 신중한 선택이 필수적이다.
너도밤나무는 전체적인 수형과 잎의 물결이 아름답고, 맹아력이 강해 해안가 경관수, 공원수 및 가로수에 뛰어난 가치를 가진다. 실제 유럽권에서는 이미 유럽너도밤나무 (Fagus sylvatica L.)의 여러 품종 개량을 통해 조경수로 활발히 이용하고 있는데, 최근 우리나라에서도 너도밤나무를 조경수로 활용하고자 하는 움직임이 보이고 있다. 다만, 지금까지 우리나라에서는 울릉도 너도밤나무 관련한 여러 연구들이 진행되어 왔으나 육지에서 너도밤나무의 생육 가능성에 대한 의문과 함께 관련 연구가 부재하였다. 이번 연구를 통해 육지에서 너도밤나무속 생육 가능성을 살펴볼 수 있었으며 너도밤나무속 3종에 대한 기후 특성을 살펴보았다. 추후 조경수 활용을 위한 기초자료가 될 것으로 판단된다.