2.1 자료 수집

본 연구에서 활용된 자료는 한강물환경연구소의 “기후변화가 수생태계에 미치는 영향과 대응전략(2010∼2012)”, 환경부의 “환경생태유량 산정기준 연구 및 시범산정(2014)”, 환경부 및 국립환경과학원의 “수생태계 건강성 조사 및 평가(2012∼2013)”, “하천 수생태계 현황 조사 및 건강성 평가(2014∼2021)”, “수생태계 참조하천 선정 및 활용방안 마련 연구(2018)”, “4대강 보 개방에 따른 수생태계 변화 조사(2018∼2020)”, 그리고 국립환경과학원의 “생물측정망 모니터링 및 평가기법 개발연구(2017)”에서 제공하는 지점 및 정점 조사 결과의 일부를 기반으로 하였다.

총 표본 단위는 23,957개, 총 조사지점은 6,664개로, 연구 대상 지역은 민통선 이남부터 제주도, 서해안 접경부터 동해안 접경, 해발고도 0∼888 m까지 전국적인 범위를 포함한다. 각 조사지점에서 수집된 자료는 하루살이속 약충의 출현 여부와 해당 지점의 물리⋅화학적 환경인자인 서식지 유형, 미소서식처 및 수질에 관한 것이며, 성충에 대한 조사는 포함되지 않았다. 또한, 일부 분석에서는 환경 요인 자료가 없는 지점을 제외하였다. 이는 생물과 환경 요인 간의 상관성을 분석하기 위한 조치로, 최종적으로 분석에 포함된 조사지점은 4,962개이다. 이들 지점에서는 연중 또는 연간 반복 조사가 이루어졌다.

하루살이속 중 사할린하루살이는 약충 단계에서 동양하루살이와 형태적으로 구분이 어렵다. 이에 따라 본 연구에서 사용된 대부분의 자료에서는 사할린하루살이가 동양하루살이에 포함되어 기록되었다. 또한, 두 종은 동소서식하는 것으로 보고된 바 있으며^{(Hwang et al., 2013)}, 이에 따라 본 연구에서는 동양하루살이-사할린하루살이 군으로 분석을 진행하였다.

2.2 현장조사

본 연구에서 고려한 물리화학적 환경 요인 중 서식지는 고도(Altitude, m), 미소서식처는 수심(Depth, m), 평균 유속(Mean current velocity, cm/s) 및 평균중위입경(Mean diameter, Φ)이 포함되었으며, 수질 요인으로는 pH, 용존산소(Dissolved Oxygen, DO, mg/L), 5일 생물학적 산소요구량(5-day Biochemical Oxygen Demand, BOD₅, mg/L), 총부유물질(Total Suspended Solids, TSS, mg/L), 총질소(Total Nitrogen, TN, mg/L), 총인(Total Phosphorus, TP, mg/L)이 포함되었다.

수심이 얕은 4,682개 지점은 쇠자를 이용하여 측정되었고 수심이 깊은 280개 지점은 센서가 내장된 휴대용 수심측정기로 측정된 자료가 수집되었다. 유속 측정은 수심이 낮은 지점에서 센서를 직접 설치하여 유속을 측정하였고, 수심이 깊은 곳에서 부유식 유속계를 활용하였다. 자유 흐름 상태에서 하천의 직선 구간과 곡류 구간, 또는 바닥과 수면 사이에 따른 차이를 고려하여 유속과 수심의 관계를 나타내는 특정 에너지(Specific Energy), 프로우드 수(Froude Number), 레이놀즈 수(Reynolds Number) 및 유체 경계층(Boundary Layer)을 복합적으로 평가하는 개념을 적용하였다^{(Craig, 1987)}. 이를 통해 산출된 값은 유속 변화를 반영하는 지표로 활용되었다.

하상기질 자료는 격자망을 이용하거나 육안으로 입도 분류 체계에 따라 분류되었으며^{(Wentworth, 1922)}, 각 구간별 입도의 비율을 활용하여 평균 중위 입경을 산정하였다. 수질 자료는 환경부 물환경측정망의 수질측정자료를 이용하였다. 저서성 대형무척추동물의 채집에는 수심이 얕은 지점에서는 Surber net 또는 D-frame net을 사용하였으며, 수심이 깊은 곳에서는 Dredge, Ekman grab, Ponar grab 등을 활용하여 시료를 확보하였다. 채집된 개체는 종 수준에서 동정한 후 개체 수를 기록하였다.

2.3 통계분석

본 연구에서 데이터의 통계 분석 및 시각화는 파이썬(Python)의 주피터 노트북(Jupyter Notebook, 7.2.2) 환경에서 수행하였다.

3.1 환경 변수 간의 피어슨 상관계수 클러스터링 결과

전체 조사 지점에서 수집된 환경 변수 간 상호관계 구조를 탐색하고자 피어슨 상관계수에 따른 계층적 클러스터링 분석을 수행한 결과(Fig. 1), 유의미한 클러스터 구조가 도출되었다. 수질변수인 TSS, BOD5, TN, TP는 같은 클러스터를 이루었으며 특히 TP와 TSS는 가장 인접하였다. 일반적으로 TP와 TSS는 강우와 깊은 관계를 가진 변수로 유역으로부터 강우유출 과정에 인이 토양입자에 흡착되어 수계로 유입되는 현상이 반영된 것으로 판단된다^{(Seo and Kim, 2016)}. 토양입자에 흡착되어 유출되는 인과 달리 질소는 토양의 공극수에 용존된 상태에서 약한 강우강도에도 쉽게 유출되는 특성이 있으므로 주요 수질변수 중에는 외군으로 나타난 것으로 보인다^{(Heathwaite et al., 2000; Mulholland, 1992)}. 수질변수 중 용존산소는 BOD5나 질소 및 인과는 달리 물리적 요인과 클러스터를 이루었는데, 이는 유기물의 분해나 질산화 및 영양물질에 의한 부착조류의 증식에 따른 영향보다는 고도에 따른 수온과 유속의 영향이 용존산소에 더욱 큰 영향을 미치고 있음을 시사한다^{(Kim and Oh, 2021; Moon and Kim, 2019)}. 하천의 pH 역시 물리적 요인과 클러스터를 이루었는데, 하류로 갈수록 유속이 완만해지고 영양물질의 존재 형태가 가용성으로 변하여 부착조류가 증가하고 이에 따른 광합성으로 pH가 상승하는 경향이 반영된 것으로 보인다(Mulholland, 1992; ^{Odum, 1956)}. DO와 pH는 이러한 하천의 연속적인 변화와 연관되어 있을 뿐만 아니라 부착조류의 광합성과 호흡 과정에 동시에 증감하는 변수이므로 인접한 클러스터를 이룬 것으로 판단된다.

평균 입경과 수심 역시 동일한 클러스터를 형성하였는데 일반적으로 수심이 얕은 곳은 조립질 하상을 이루는 반면, 수심이 깊은 곳은 세립질 하상이 주를 이루는 경향을 반영하고 있는 것으로 보인다^{(Allan and Castillo, 2007)}. 이러한 물리적 요인은 고도와 유속 또는 TP 등 수질요인과 다른 클러스터를 이루었는데, 고도가 낮은 저지대에도 얕은 수심의 소하천이 많이 분포하기 때문인 것으로 판단된다^{(Montgomery and Buffington, 1997; Vannote et al., 1980)}. 한편, 유속과 고도는 동일한 군집으로 분류되었으며, 이는 하천 상류로 갈수록 고도가 높아지고 유속이 증가하는 전형적인 지형-수리학적 관계를 반영하는 것으로 해석된다. 실제로 ^{Kong and Song (2023)}은 확률모형 기반 분석을 통해 동양-사할린하루살이는 느린 유속(rheophobic), 가는무늬하루살이는 빠른 유속(stenomesorheophilic)에 뚜렷한 출현 특성을 보인다고 보고된 것을 참고하여, 본 연구의 조사지점은 여울(riffle), 흐름(run), 소(pool) 등을 포함한 다양한 미소서식처 유형을 포함하였으며, 유속 측정은 하루살이 개체의 실제 출현 지점에서 수행된 값을 기반으로 수행하였다. 따라서 본 연구에서 분석한 유속 변수는 일반적인 하천 평균 유속이 아닌, 개체가 서식하는 미소 환경 특성을 반영한 값으로 해석할 수 있다.

Fig. 1. Hierarchical clustering dendrogram of environmental variables based on Pearson's correlation coefficients.

3.2 BIOENV 분석 결과

BIOENV 분석을 통해 하천 내 생물 종 조성에 영향을 미치는 환경 변수를 평가한 결과, 고도가 가장 높은 양의 상관관계(Pearson’s ρ, 0.1712)를 나타내어, 군집 구조에 가장 큰 영향을 미치는 주요 환경 변수로 확인되었다(Table 1). 그 외 입경(0.1057), pH(0.0575), 평균 유속(0.0556), DO(0.0474), 수심(0.0461), TN(0.0430), TSS(0.0381), BOD5(0.0093) 등의 변수는 양의 상관관계를 나타냈으나, 그 크기는 비교적 낮아 생물 종 분포에 대한 영향은 상대적으로 제한적인 것으로 해석된다. 한편, TP은 유일하게 Pearson 상관계수 -0.0039로 음의 상관관계를 나타냈으며, 생물 종 구성과의 직접적인 연관성은 미약한 것으로 분석되었다. 전체적으로 상관계수 값이 0.2 이하로 나타난 점은, 단일 환경 변수보다는 다수의 변수들이 복합적으로 작용하여 생물 군집 구조를 형성하고 있음을 시사한다. 이는 개별 변수의 독립적인 기여도 외에도 복수의 환경 요인이 함께 작용하여 생물 군집 구조를 결정한다는 점을 강조하며, 세 종의 분포 차이를 결정하는 환경 반응성의 구조적 맥락을 포괄적으로 분석했다는 점에서 의의가 있다.

결론적으로, 하루살이속 종들은 환경변수의 절대값보다는 범위 및 조합, 그리고 군집구조 등의 복합적인 요인에 의해 출현이 조절되며, 종마다 수질 인자에 대한 내성 수준 및 미소 서식처의 차이에 따라 상대적인 출현 빈도가 달라질 수 있기 때문에, 오염 관련 변수(TN, TSS, BOD5 등)의 공동 작용 가능성과 하천 구조 및 지형 요인(고도, 유속, 입경, 수심 등)의 물리적 맥락을 해석에 포함하는 것이 필요하다. 단일 변수 중심의 BIOENV 결과를 생태적으로 해석하기 위해서는 이러한 다층적 특성과 종간 상호작용을 고려한 통합적 분석이 필요하며, 개별 종 또는 종군 단위로의 분리 분석과 더불어 각 종의 생태적 특성과 내성 범위를 고려한 다층적 접근이 요구된다.

3.3 중복분석 결과

중복분석 결과, 고도와 유속은 두 축 모두에서 주요한 환경 구배 축을 형성하며, 하루살이속 종 들의 공간적 생태 지위를 구분하는 핵심 요인으로 작용하는 것으로 나타났다(Fig. 2). 특히 중복분석 기반 환경 설명력(R²)은 무늬하루살이(0.1227)가 가장 높게 나타났으며, 동양하루살이-사할린하루살이(0.0421), 가는무늬하루살이(0.0191) 순으로 점차 감소하였다. 이는 종별로 환경 구배에 대한 민감도 및 생태적 전략에 차이가 존재함을 시사한다. 중복분석 축의 분산 설명력은 Axis 1이 66.94 %, Axis 2가 32.25 %로, 두 축이 전체 종 분포의 99.19 %를 설명하였다. 특히 Axis 1은 고도 및 유속을 중심으로 한 청정 환경 구배를, Axis 2는 TSS, BOD5, TP 등의 오염 및 서식처 구조 관련 변수 구배를 반영하는 것으로 해석된다. 이 두 차원을 통해 하루살이속 종들이 환경 축에 따라 뚜렷이 구분되는 생태적 지위를 나타냈다.

동양하루살이-사할린하루살이는 TP, BOD5, TSS, 수심, 입경 등 오염 및 미소 서식처 관련 지표들과 유사한 방향성을 보이며, 고도 벡터와는 반대 방향에 위치하였다(Fig. 2a). 이러한 벡터 배열은 해당 종이 다양한 오염도 및 서식 환경에 대해 내성을 가지는 광서식성(eurytopic) 종임을 시사하며, 실제 출현 점들도 중복분석 공간 전반에 고르게 분포하여 이러한 특성과 일치한다^{(Park and Park, 1999)}. 수질 환경의 다양성에 대응 가능한 생리적 유연성을 지닌 종으로 해석되며, 이는 내오염성 지표종 또는 참조 지점 설정 시 보완지표로 활용 가능함을 의미한다. 또한, 이 종의 중심 좌표는 (0.0033, –0.0016)으로 중복분석 공간의 중심에 가장 가까이 위치하여, 특정 환경 구배보다는 다양한 조건에 대한 분포 가능성을 시사한다.

반면, 무늬하루살이는 고도 및 유속 벡터에 평행하게 정렬되며, TP, TSS, BOD5 등 오염 지표들과는 반대 방향에 위치하였다(Fig. 2b). 높은 환경 설명력(R², 0.1227)과 분포의 공간 집중도는 이 종이 청정한 고지대의 빠른 유속 환경에 특화된 민감성 종임을 시사한다. 실제로 Axis 1에서의 강한 양의 기여는 이 종이 상류 청정 수역을 대표하는 생태 지표종으로 기능할 수 있음을 보여준다. 중심 좌표는 (–0.0470, 0.0225)로, 고도 및 유속 구배 축의 양의 방향에 강하게 정렬되어 있음이 수치적으로도 확인된다.

가는무늬하루살이는 고도 및 유속 벡터와 유사한 방향에 위치하나, 중심축에 가까운 분포를 보이며, 벡터의 길이 또한 짧고 수렴적인 특성을 보였다(Fig. 2c). 이는 좁은 환경 구간에 국한된 협서식성(stenotopic) 종으로 해석되며, 환경 구배에 대한 제한적 반응성을 나타낸다. 수질 오염 지표들과는 명확히 분리되어 있으며, 분포 범위가 협소하고, 특정 미세 서식처(예: 얕고 고요한 고지대 흐름)에서만 나타나는 특성이 강하게 나타났다. 중심 좌표 또한 (–0.0415, 0.0173)으로 무늬하루살이와 유사하되, 분포의 밀도는 낮고 집중도는 떨어지는 경향을 보인다.

이와 같이, 중복 분석 결과는 하루살이속 종들이 동일한 분류군 내에서도 환경 구배에 따라 뚜렷한 생태적 반응성과 공간적 분화를 보인다는 점을 시사한다. 특히 종별 중심 좌표와 축별 설명력, 환경 벡터의 방향성과 일치하는 공간 분포 특성은 단순한 시각적 해석을 넘어 수치 기반의 생태적 차별성을 뒷받침하며, 생물지표 개발이나 환경 평가 기준 설정 시 종 단위의 접근이 중요함을 강조한다.

Fig. 2. Redundancy analysis (RDA) ordination plots showing relationships between environmental variables and each Ephemera species: (a) E. orientalis–sachalinensis, (b) E. strigata, and (c) E. separigata.

3.4 다변량 회귀 모델 결과

부분최소제곱회귀(PLSR)와 일반최소제곱회귀(OLS) 모델을 적용하여 비교하였다(Fig. 3). 두 모델의 성능 평가 결과, RMSE 값은 PLSR : 0.1146, OLS : 0.1154로 나타나 두 모델 간 성능 차이는 미미했으나, 환경 변수에 대한 공통적 경향성을 보였다는 점에서 결과의 신뢰성을 높일 수 있었다.

두 모델 모두에서 고도는 상위 수준의 영향력을 가지는 주요 변수로 도출되었으며, 고도에 따른 하루살이 종의 출현이 뚜렷하게 구분되는 경향을 보였다. 동양하루살이-사할린하루살이는 고도와 음의 관계를 보였으며, PLSR : -0.0086, OLS : -0.0162로 고도가 낮은 지역에 더 많이 출현하는 경향을 보였다. 반면, 무늬하루살이 PLSR : 0.0074, OLS : 0.0125, 가는무늬하루살이 PLSR : 0.0015, OLS : 0.0040으로 고도가 높을수록 출현율이 증가하는 경향을 보였다. 실제, 고도에 따른 출현 특성을 확률분포모형으로 분석한 연구^{(Kong and Kang, 2023)}에 따르면, 동양하루살이-사할린하루살이 군은 저고도(평균 251 m, 최빈값 124 m)에 폭넓게 분포하며 열선호성(thermophilic)을 지닌 반면, 무늬하루살이는 중위고도(평균 474 m, 최빈값 492 m)의 중간적 온도대에서 출현하는 중온성(mesophilic), 가는무늬하루살이는 높은 고도(평균 620 m, 최빈값 760 m)의 냉수성(thermophobic) 서식처에 특화된 양상을 보인다 .

유속은 무늬하루살이 PLSR : 0.0049, OLS : 0.0064로 양의 상관성을 보여, 고도가 높은 지역의 빠른 유속을 선호하는 경향을 보였다. 동양하루살이-사할린하루살이(PLSR : 0.0071, OLS : 0.0029)는 양의 값을 보여 약한 상관을 나타냈으며, 가는무늬하루살이(PLSR : 0.0001, OLS : -0.0019)는 일관된 경향을 보이지 않았다. 수심은 고도와 반대되는 패턴을 보였으며, 동양하루살이-사할린하루살이는 PLSR : 0.0147, OLS : 0.0171로 양의 상관을 나타냈다. 이는 깊은 수심에서 높은 출현율을 보이는 경향과 일치하며, 무늬하루살이(PLSR : -0.0039, OLS : -0.0044) 및 가는무늬하루살이(PLSR : -0.0015, OLS : 0.0002)는 수심에 대해 음 또는 미약한 반응성을 보였다.

유속⋅수심⋅입도에 기반한 물리적 미소서식처 적합도 연구^{(Kong and Song, 2023)}에서는 동양하루살이-사할린하루살이는 깊은 수심, 느린 유속, 세립질 하상에 적응된 광호심성(eurybathophilic) 및 광쇄석성(eurypsephophilic) 특성을 가지며, 무늬하루살이는 얕은 수심, 빠른 유속, 조립질 하상을 선호하는 광유수성(euryrheophilic), 암석성(lithophilic) 종으로 분류된다. 가는무늬하루살이는 낮은 수온과 깨끗한 수질에 특화된 청정 협서식성을 보여 본 연구와도 일치 하였다.

수질 오염 관련 변수(TSS, TP, BOD5)는 종별로 상이한 영향을 미쳤다. 동양하루살이-사할린하루살이는 TSS(PLSR : 0.0063, OLS : 0.0089), TN(PLSR : 0.0087, OLS : 0.0075), TP(PLSR : 0.0089, OLS : 0.0209) 등과 모두 양의 상관을 보여 오염도가 높은 환경에 대한 내성을 시사하였다. 반면, 무늬하루살이와 가는무늬하루살이는 대부분 음의 상관을 보이며, 상대적으로 청정 환경에 특화된 생태 반응성을 나타냈다.

동양하루살이-사할린하루살이 군이 높은 농도의 BOD5에 적응한 광부수성(eurysaprobic) 종으로 평가되며, 높은 TP 및 TSS 농도 환경에서는 제한적인 적응성을 보이는 것으로 알려져 있다^{(Park and Park, 1999)}. 또한, ^{Jung and Kong (2023)}의 수질적합도 평가 결과에선 동양하루살이-사할린하루살이는 BOD5(0.4∼1.9 mg/L), TSS(2.5∼22.4 mg/L), TP(0.024∼ 0.086 mg/L) 등 오염도가 높은 환경에서 서식 적합도를 보였으며, 이는 본 연구에서 분석된 오염 지표들과의 정합성과 일치하였다. 반면, 무늬하루살이는 낮은 BOD5, TP, TSS 농도에서 주로 출현하는 빈부수성(oligosaprobic) 종이며, 가는무늬하루살이는 무늬하루살이와 일부 서식지를 공유하지만, 특히 낮은 BOD5 농도에서 더욱 제한적으로 분포하는 협빈부수성(stenooligosaprobic) 종으로 보고 되었다.

해외 연구 사례에서도 하루살이속 종들의 서식 분포 특성 결과는, 고도, 유속, 하상 경사 등의 환경 요인에 따라 뚜렷이 구분되며, 이러한 특정 환경 조건에 선택적으로 적응한 지위 분화(niche differentiation)를 반영한다는 결과가 보고되어^{(Okamoto et al., 2022; Okamoto and Tojo, 2021)}, 본 연구 결과와 높은 일치성을 보였다.

Fig. 3. Comparison of environmental factor effects on Ephemera species based on PLSR (a) and OLS (b) model results.

3.5 고도에 따른 상대적 출현 비율 변화

앞선 다변량 분석을 통해 하루살이속 종들의 출현과 분포에 가장 큰 영향을 미치는 주요 환경 요인이 고도임을 확인하였으며, 이를 정량적으로 검증하기 위해 실시한 일원분산분석(ANOVA) 결과, 고도 구간 간 출현 비율에서 통계적으로 유의한 차이를 보였다(p < 0.001). 이는 고도가 이들의 출현 양상에 유의미한 영향을 미친다는 점을 통계적으로 뒷받침하는 결과이다.

그러나 전체 표본 수는 저고도(0 m)에 집중되어 있었으며, 특히 0∼100 m 구간에서만 3,800개의 표본 수가 관찰되었고, 300 m 이상의 고도에서는 표본 수가 급격히 감소하는 불균형을 보였다. 이에 따라 출현 경향을 고도에 따라 공정하게 비교하기 위해 상대적 출현 비율(%)에 기반한 분석을 추가로 수행하였다(Fig. 4). 이 과정에서는 고도를 100 m 단위로 구분하여 각 구간별 평균 상대 출현 비율을 산출하였으며, 이를 통해 종 간 절대 개체 수 차이를 보정하고, 고도 변화에 따른 상대적 종 조성 변화를 보다 명확히 비교하고자 하였다.

분석 결과 동양하루살이-사할린하루살이는 저고도 지역에서 압도적으로 높은 출현 비율을 나타냈다. 0∼100m에서 96.84 %, 100∼200 m에서 81.59 %, 200∼300m에서 2.02 %의 비율을 보이며 분포의 중심이 저고도에 집중되어 있었다. 이후 고도가 증가할수록 출현 비율은 꾸준히 감소하였고, 300∼400 m에서는 67.31 %, 400∼500 m에서는 52.47 %로 감소하였다. 이후 고도가 상승할수록 급격한 감소세를 보이며 500∼600 m에서 14.36 %, 600∼700 m에서 7.17 %, 700∼900 m 고도에서는 완전히 출현하지 않았다(0.00 %).

무늬하루살이는 저고도에서 낮고, 중위 고도에서 매우 높은 출현율을 보였다. 0∼300 m 구간에서는 각각 2.62 %, 16.19 %, 15.08 %로 비교적 낮았지만, 고도가 높아질수록 빠르게 증가하여 300∼400 m에서는 29.34 %, 400∼500 m에서는 38.85 %, 500∼600 m에서는 85.64 %, 그리고 600∼700 m에서는 92.83 %로 최고치를 기록하였다. 이후 고도가 더 높아지면 출현 비율은 감소하여 700∼800 m에서 58.19 %, 800∼900 m에서는 21.24 %를 보였다.

가는무늬하루살이는 고도에 따라 점진적으로 출현율이 증가하는 분포 패턴을 나타냈다. 0∼300 m에서는 각각 0.54 %, 2.22 %, 2.90 %로 매우 낮았고, 300∼400 m에서 3.35 %, 400∼500 m에서는 8.69 %로 소폭 증가하였다. 이후 500∼700 m 구간에서는 출현이 관찰되지 않았으나, 700∼800 m에서는 41.81 %, 800∼900 m에서는 78.76 %로 급격히 증가하며 고도에 따른 뚜렷한 적응 특성을 보였다. 이는 가는무늬하루살이가 높은 고도, 빠른 유속, 낮은 수온, 양호한 수질 등의 청정 환경 조건에 특화된 생태적 특성을 뒷받침하는 결과이다.

이는 하천의 고도 변화가 종의 출현과 다양성 유지에 직접적으로 연관될 수 있음을 의미하며, 물리적 서식처 특성과 연계된 환경 요인이 생물다양성 유지에 중요한 역할을 한다는 기존 연구 결과와 일치한다(Kong and Kang, 2023; ^{Lee et al., 1995; Okamoto et al., 2022; Okamoto and Tojo, 2021; Park et al., 1997)}.

Fig. 4. Relative abundance (%) of Ephemera species across altitude groups (100 m intervals), based on averaged normalized data.