Radulesco Thomas
email: thomas.radulesco@ap-hm.fr

Lionel Meister

Gilles Bouchet

Jérome Giordano

Patrick Dessi

Pierre Perrier

Justin Michel

Meister Lionel

Bouchet Gilles

Giordano Jérome

Dessi Patrick

Perrier Pierre

Michel Justin

Thomas Radulesco

Functional relevance of computational fluid dynamics in the field of nasal obstruction: A literature review

Keywords: CFD, nasal airway, obstruction, surgery, functional, septoplasty, rhinoplasty, empty nose syndrome, nasal valve, turbinectomy

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

Nasal airway obstruction (NAO) is a common symptom affecting the quality of life of patients [START_REF] Rhee | Quality of life assessment in nasal airway obstruction[END_REF]. It may involve different etiologies including septal deviations (SD), nasal valve collapse (NVC), or inferior turbinal hypertrophy (ITH) [START_REF] Do | Predicting postsurgery nasal physiology with computational modeling: current challenges and limitations[END_REF]. The clinical evaluation of NAO can be subjective, by self-questionnaires (NOSE, SNOT 22) and Visual Analog Scale (VAS), or reported by different tests such as rhinomanometry or Peak Nasal Inspiratory Flow (PNIF) (3) (4) [START_REF] Pendolino | Comparison between unilateral PNIF and rhinomanometry in the evaluation of nasal cycle[END_REF].

Computational fluid dynamics (CFD) consists in studying fluid movements, or its effect, by numerical resolution of equations governing this fluid. It can be applied to nasal anatomy: over the past 10 years, there have been many studies focusing on CFD in nasal airway [START_REF] Burgos | A CFD approach to understand nasoseptal perforations[END_REF] [START_REF] Chen | Numerical simulation of the effects of inferior turbinate surgery on nasal airway heating capacity[END_REF]. CFD can accurately model the airflow and its conditioning in the nasal cavities and calculate new variables, unknown in the world of medicine, such as Heat Flux (HF) or Wall Shear Stress (WSS) [START_REF] Dayal | Impact of Middle versus Inferior Total Turbinectomy on Nasal Aerodynamics[END_REF]. 3D reconstruction tools can help the clinician better understand nasal anatomy, and CFD data to better analyze NAO genesis [START_REF] Kimbell | Changes in nasal airflow and heat transfer correlate with symptom improvement after surgery for nasal obstruction[END_REF]. The final goal is to guide surgical procedures. However, these data are not often correlated with patients perceptions or measurements carried out in clinical routine [START_REF] Rhee | Role of virtual surgery in preoperative planning: assessing the individual components of functional nasal airway surgery[END_REF]. The correlation between clinical evaluation of nasal breathing and CFD is essential to improve the interpretation of the latter. There is a lack in literature of comparative study analyzing the different results obtained, as Leite et al. reports in his recent literature review. [START_REF] Otto | Computational fluid dynamics evaluation of posterior septectomy as a viable treatment option for large septal perforations[END_REF] [START_REF] Leite | The clinical implications of computerised fluid dynamic modelling in rhinology[END_REF] We conducted a literature review of CFD studies regarding the nasal airway, focusing on articles comparing clinical data to CFD. The goal is to provide what are the most reliable CFD measures: this study investigates the correlation between CFD and NAO.

METHODS

The selected studies were obtained from the database of the US National Library of Medicine (PubMed) online database, MEDLINE (Ovid), Google scholar, and Cochrane Library using a combination of MeSH terms (nose, paranasal sinus, fluid dynamics, rhinology) and no MeSH terms (CFD, nasal airway, nasal airflow, numerical, nasal symtoms). Studies that do not incorporate objective or subjective clinical assessment were excluded.

RESULTS

Articles selection

Table 1 reports all selected studies and the ways they assess clinical evaluation. Of the 258 articles selected, 47 dealt with nasal airway and corresponded to our topic. Thirteen of these 47 articles were found to be revevant, including with clinical measurements (Figure 1). This review sums up the correlation between clinical assessment of NAO and CFD results. We did not apply a time-out filter for the research.

Nasal Airway obstruction: Clinical Assessment Clinical assessment of NAO can be performed by many tools. Regarding subjective assessment, most authors use the NOSE questionnaire and / or the VAS scale [START_REF] Rhee | Role of virtual surgery in preoperative planning: assessing the individual components of functional nasal airway surgery[END_REF] [START_REF] Floyd | Systematic Review and Meta-analysis of Studies Evaluating Functional Rhinoplasty Outcomes with the NOSE Score[END_REF]. NOSE questionnaire is a disease-specific quality-of-life instrument for NAO that has been validated in literature [START_REF] Hsu | Evaluation of nasal patency by visual analogue scale/nasal obstruction symptom evaluation questionnaires and anterior active rhinomanometry after septoplasty: a retrospective one-year follow-up cohort study[END_REF]. NOSE is a graduated 20-point scale, which may be multiplied by 5 to give a final total out of 100. A score of 100 indicates complete nasal obstruction. A NOSE score under 25 is considered normal; a score between 25 and 50 denotes low nasal obstruction and a score above 50 severe nasal obstruction. Others use the SNOT 22 questionnaire, or more specific questionnaires as in the Empty Nose Syndrome (ENS) with the ENS6Q questionnaire [START_REF] Stewart | Development and validation of the Nasal Obstruction Symptom Evaluation (NOSE) scale. Otolaryngol--Head Neck[END_REF] [START_REF] Erskine | SNOT-22 in a control population[END_REF].

Regarding the paraclinical evaluation, the most used tool remains the rhinomanometry which measures pressure and flow during normal inspiration and expiration through the nose. (4).

Acoustic rhinomatry is a measurement of cross sectional area and length of the nose and the nasal cavity through acoustic reflections. It can be correlated to minimal cross section area (MCAs) and is also used by some authors for functional assessment, although it is actually a morphological examination [START_REF] Velasquez | The Empty Nose Syndrome 6-Item Questionnaire (ENS6Q): a validated 6-item questionnaire as a diagnostic aid for empty nose syndrome patients[END_REF]. We did not found any study involving PNIF.

CFD Calculations

The main data measured by CFD are as follows:

-Airflow (AF) (uni-or bilateral, in mL/s), and Airflow partitioning (%) corresponding to the balance between obstructed and non-obstructed nasal fossa.

-Nasal resistance (CFD-NR) to uni-or bilateral inspiration and expiration, (Pa/(mL/s)):

Nasal resistance (ΔP/QV) is defined as the ratio of the drop in transnasal pressure ΔP (nostrils to nasopharynx or to the end of the septum) to the volumetric airflow rate.

-Total Heat flux (HF) (W/m2), which is the rate of heat transfer across a surface per unit time and area and measures heat loss from the nasal mucosa to the inspired air; Surface Area where Heat Flux > 50 W/m 2 (SAHF50); or Peak Heat Flux defined as the value above which only 1 cm2 of mucosa is exposed to.

-Total Pressure (Pa) and pressure drop -Air velocity (m/s) and Maximal velocity (Vmax) -Wall shear stress (WSS) (Pa), which is a friction force generated when moving air contacts the nasal walls -Streamlines, allowing airway visualization The decision to analyze all or part of these data depends on the purpose and design of the study. Some data may be concerned only with airflow (velocity, pressure), and some only with air conditioning.

Correlation between CFD and clinical assessment

Airflow (AF)

Casey et al. Found a strong correlation between unilateral nasal AF and NOSE (r=-0.55, p=0.0016) or VAS (r=-0.45, p=0.0056) [START_REF] Geißler | Rhinomanometry and acoustic rhinometry[END_REF]. It is the same for Gaberino et al., but only after correction of the nasal cycle: without correction of the latter, AF was correlated with NOSE (p = 0.048), but not with VAS [START_REF]Casey_Correlation between Subjective Nasal Patency and Intranasal Airflow Distribution[END_REF]. Kimbell also reports a moderate correlation between NOSE and AF [START_REF] Gaberino | Estimates of nasal airflow at the nasal cycle mid-point improve the correlation between objective and subjective measures of nasal patency[END_REF].

In patients improved after surgery for NAO, several studies found improvement of AF in the narrow side, and a deterioration of AF in the non-narrow side after surgery [START_REF] Gaberino | Estimates of nasal airflow at the nasal cycle mid-point improve the correlation between objective and subjective measures of nasal patency[END_REF][START_REF] Kimbell | Changes in nasal airflow and heat transfer correlate with symptom improvement after surgery for nasal obstruction[END_REF][START_REF] Kim | Correlation between nasal airflow characteristics and clinical relevance of nasal septal deviation to nasal airway obstruction[END_REF].

CFD-calculated Nasal Resistances (CFD-NR)

Nasal resistances are among the most analyzed data in the literature concerning NAO. By analyzing two groups of patients (NAO and non-NAO), Casey et al. found a statistically significant difference (p = 0.0006) in CFD-NR measurements [START_REF] Geißler | Rhinomanometry and acoustic rhinometry[END_REF]. However, in his study on 30 patients, there was no correlation between CFD-NR and the NOSE or VAS. It is the same for Gaberino et al. who need to correct the nasal cycle to obtain a correlation between CFD-NR with the clinical data. Thus, and only after correction of the nasal cycle, CFD-NR were correlated with NOSE and VAS (r = 0.55, p = 0.005 and r = -0.58, p = 0.003 respectively) [START_REF]Casey_Correlation between Subjective Nasal Patency and Intranasal Airflow Distribution[END_REF]. Kimbell found that CFD-NR were moderately correlated with NOSE and strongly with VAS, similar results for Rhee et al. [START_REF] Gaberino | Estimates of nasal airflow at the nasal cycle mid-point improve the correlation between objective and subjective measures of nasal patency[END_REF][START_REF] Kim | Correlation between nasal airflow characteristics and clinical relevance of nasal septal deviation to nasal airway obstruction[END_REF][START_REF] Kimbell | Computed nasal resistance compared with patient-reported symptoms in surgically treated nasal airway passages: a preliminary report[END_REF]. In contrast, Zhao et al. found no correlation between CFD-NR and VAS [START_REF] Rhee | Role of Virtual Surgery in Preoperative Planning: Assessing the Individual Components of Functional Nasal Airway Surgery[END_REF]. Concerning the non-narrow side, it is possible to find an increase in postoperative resistance [START_REF] Zhao | Regional peak mucosal cooling predicts the perception of nasal patency[END_REF].

Regarding objective measurements, nasal resistances obtained with rhinomanometry and CDF-NR are correlated with no significative differences (p>0.05 and r=0.41, p<0.05) [START_REF] Rhee | Role of Virtual Surgery in Preoperative Planning: Assessing the Individual Components of Functional Nasal Airway Surgery[END_REF][START_REF] Hildebrandt | Evaluation of the intranasal flow field through computational fluid dynamics[END_REF].

Total Heat Flux (HF), SAHF50, Peak Heat Flux

HF analysis have now a major place in CFD in the field of NAO. NOSE et VAS were strongly correlated with total unilateral HF (r=-0.48, p= 0.0075 and r=0.43, p=0.0166 respectively), and with SAHF50 (r=-0.55, p=0.0016 and r=-0.51, p=0.0038 respectively) in Casey's study. The correlation between HF and NOSE was moderate [START_REF] Geißler | Rhinomanometry and acoustic rhinometry[END_REF]. Sullivan et al. also found that HF and SAHF50 was strongly correlated with NOSE (r=0.76, p<0.01) and VAS (r=0.5, p<0.05) (28).

Gaberino et al. found a statistical significative difference between pre and post-operative

analyze for HF (p=0.027); after correction of the nasal cycle, NOSE was correlated to HF and SAHF50 but VAS was only correlated to HF [START_REF]Casey_Correlation between Subjective Nasal Patency and Intranasal Airflow Distribution[END_REF].

For Zhao et al., only the peak nasal mucosal heat loss posterior to the nasal valve was correlated to clinical impairment (r=-0.46, p<0.01) [START_REF] Rhee | Role of Virtual Surgery in Preoperative Planning: Assessing the Individual Components of Functional Nasal Airway Surgery[END_REF].

Pressure

Hildebrandt et al. shows a diminution of pressure drop in the nasal valve area after surgery for NAO (Figure 3a) [START_REF] Zhao | Regional peak mucosal cooling predicts the perception of nasal patency[END_REF]. It highlights that the nasal valve is a key-zone for NAO. When comparing 3 symptomatic versus asymptomatic patients, Kim et al. also found a pressure drop more important in patients with NAO (p<0.05) [START_REF] Kimbell | Changes in nasal airflow and heat transfer correlate with symptom improvement after surgery for nasal obstruction[END_REF]. (Figure 3b) Maximum Velocity (Vmax) For Kim et al. Vmax was higher in patients with SD (p=0.05) [START_REF] Kimbell | Changes in nasal airflow and heat transfer correlate with symptom improvement after surgery for nasal obstruction[END_REF]. Vmax could also be correlated to clinical impairment: when the SD is caudal, Vmax and symptoms are higher than in posterior SD [START_REF] Sullivan | Perception of better nasal patency correlates with increased mucosal cooling after surgery for nasal obstruction[END_REF].

Wall Shear Stress (WSS)

WSS in nasal valve area was higher in patients with SD (p=0.05) [START_REF] Kimbell | Changes in nasal airflow and heat transfer correlate with symptom improvement after surgery for nasal obstruction[END_REF]. In patients with ENS, WSS is lower in the inferior region, due to absence of inferior turbinate [START_REF] Liu | Effects of septal deviation on the airflow characteristics: using computational fluid dynamics models[END_REF][START_REF] Li | Computational fluid dynamics and trigeminal sensory examinations of empty nose syndrome patients[END_REF].

Streamlines

Regarding nasal airway, patients with NAO have less air in the middle pathway (31 +-18 mL/s vs 68 +-10 mL/s, p<0.001). It was strongly correlated with NOSE (r=-0.76, p<0.001) and VAS scores (r=0.9, p<0.0001) [START_REF] Geißler | Rhinomanometry and acoustic rhinometry[END_REF].

In a more specific population with ENS patients, there's a surprising modification in airflow when inferior turbinate is missing : we observe a frank diminution of airflow in the inferior region and an augmentation of the latter in middle region (p<0.05) [START_REF] Li | Computational fluid dynamics and trigeminal sensory examinations of empty nose syndrome patients[END_REF].

Nasal geometry

No differences were obtained by acoustic rhinometry and by 3D measurement (all p > 0.05) [START_REF] Hildebrandt | Evaluation of the intranasal flow field through computational fluid dynamics[END_REF].

DISCUSSION

Synopsis of new findings

One of the difficulties regarding outcomes of esthetic and/or functional surgery is the objectivity of the evaluations carried out and their correlation with the patient's perception. [START_REF] Li | Investigation of the abnormal nasal aerodynamics and trigeminal functions among empty nose syndrome patients[END_REF]. While self-evaluation using quality of life questionnaires remains the gold standard, results of objective tests such as rhinomanometry remain controversial [START_REF] Barone | A systematic review of patient-reported outcome measures after rhinoplasty[END_REF]. The development and improvement of CFD processes should allow for reliable and correlated assessments of patientreported outcomes. The advantage of CFD lies not only in airflow analysis but also in its conditioning: e.g. in heat flux measurements that could be influental in the etiopathogeny of NAO [START_REF] Pawar | Objective measures in aesthetic and functional nasal surgery: perspectives on nasal form and function[END_REF]. This is very innovative as these data are not measurable by conventional investigations.

Functional relevance of CFD

The study of CDF-NR and AF is easy because these physical data have been used for a long time in rhinology. These are values that are commonly used by clinicians.

Regarding the comparison between patient perceptions and nasal resistances, authors agree that unilateral evaluation is more reliable than bilateral evaluation [START_REF] Kim | Correlation between nasal airflow characteristics and clinical relevance of nasal septal deviation to nasal airway obstruction[END_REF] [START_REF] Zhao | Perceiving nasal patency through mucosal cooling rather than air temperature or nasal resistance[END_REF] [START_REF] Naito | The effect of L-menthol stimulation of the major palatine nerve on subjective and objective nasal patency[END_REF]. Thus, it is necessary to focus on narrow-side and non narrow-side rather than total AF or total resistance [START_REF]Casey_Correlation between Subjective Nasal Patency and Intranasal Airflow Distribution[END_REF]. While in most studies there is a significant improvement after surgery on the narrow-side, contralateral deterioration is often found, probably because of medialization of the nasal septum. However, these patients do not have a functional complaint [START_REF] Zhao | Regional peak mucosal cooling predicts the perception of nasal patency[END_REF]. It therefore appears that the balance of airflow partitioning is important to consider in the functional complaint.

Airflow partitioning is easily evaluable with CFD.

In the field of nasal resistances, we find globally similar results between clinical nasal resistances and CFD-NR. Some studies have already shown the correlation between clinical measures of nasal resistances and CFD-NR (27) (25) [START_REF] Naito | The effect of L-menthol stimulation of the major palatine nerve on subjective and objective nasal patency[END_REF]. However, these studies focus on a small number of subjects or on a healthy population. More studies need to be done about it.

The study of correlations between CFD-NR with patient perception finds heterogeneous results. Kimbell reports a moderate correlation between CFD-NR and NOSE or VAS [START_REF] Kimbell | Changes in nasal airflow and heat transfer correlate with symptom improvement after surgery for nasal obstruction[END_REF][START_REF] Kimbell | Computed nasal resistance compared with patient-reported symptoms in surgically treated nasal airway passages: a preliminary report[END_REF]. Gaberino et al. found that only AF and NOSE are correlated before correction of the nasal cycle in their population, results similar in other studies [START_REF]Casey_Correlation between Subjective Nasal Patency and Intranasal Airflow Distribution[END_REF] [START_REF] Geißler | Rhinomanometry and acoustic rhinometry[END_REF][START_REF] Rhee | Role of Virtual Surgery in Preoperative Planning: Assessing the Individual Components of Functional Nasal Airway Surgery[END_REF]. Many reasons can explain these discordant results. Firstly, the sensation of NAO is not only due to the increase of nasal resistances. This is why clinical measurements such as rhinomanometry are themselves discussed. Indeed, several studies have shown that rhinomanometry is not always correlated with clinical impairment, and more complex factors intervene in the perception of NAO. [START_REF] Zachow | Visual Exploration of Nasal Airflow[END_REF] [START_REF] Chandra | Diagnosis of nasal airway obstruction[END_REF]. Secondly, it is necessary to pay attention to the nasal cycle in CFD study, which can be a confounding factor [START_REF] Umihanic | The Discrepancy Between Subjective and Objective Findings After Septoplasty[END_REF] [START_REF] Kahana-Zweig | Measuring and Characterizing the Human Nasal Cycle[END_REF] [START_REF]Casey_Correlation between Subjective Nasal Patency and Intranasal Airflow Distribution[END_REF]. If a significant difference in the volume of the inferior turbinates is noted on CTscans, it is necessary to make a correction thanks to virtual surgery, or to take it into account in the interpretation of the results. [START_REF]Casey_Correlation between Subjective Nasal Patency and Intranasal Airflow Distribution[END_REF]. Third, tissue compliance is currently not considered in CFD studies, working on rigid walls. It is recognized that NAO problems due to nasal valve may be static or dynamic [START_REF] Patel | Simulating the Nasal Cycle with Computational Fluid Dynamics[END_REF] [START_REF] Bloching | Disorders of the nasal valve area[END_REF]. Dynamic problems can be highlighted by an examination called rhinoresisometry [START_REF] Aksoy | Role of nasal muscles in nasal valve collapse[END_REF]. If static (architectural) problems are well explored in CFD, simulating a NVC appearing during breathing would require modeling deformable walls according to the airflow. We did not find major difference between NOSE and VAS, their degree of correlation varying according to the studies. We recommend to use both, as well as more specific questionnaires like the ENS6Q for the ENS. In addition, we did not find any CFD study with PNIF as an objective measurement of NAO.

Recent studies show the major interest that seems to represent HF in rhinology. [START_REF] Geißler | Rhinomanometry and acoustic rhinometry[END_REF][START_REF]Casey_Correlation between Subjective Nasal Patency and Intranasal Airflow Distribution[END_REF][START_REF] Liu | Effects of septal deviation on the airflow characteristics: using computational fluid dynamics models[END_REF][START_REF] Li | Computational fluid dynamics and trigeminal sensory examinations of empty nose syndrome patients[END_REF][START_REF] Temmel | Rhinoresistometry versus rhinomanometry--an evaluation[END_REF]. For a long time, a link has been found between perception of cold and NAO perception [START_REF] Inthavong | Geometry and airflow dynamics analysis in the nasal cavity during inhalation[END_REF]. The new CFD tools allow a fine analysis of this thermal condition. In the studies reviewed, all data agree to give a major place to HF in the genesis of NAO perception. Different versions of HF can be used (Total HF, SAHF50, Peak HF), each of which has shown a statistically significant correlation with patients' perceptions (Figure 2) [START_REF] Lu | Accuracy evaluation of a numerical simulation model of nasal airflow[END_REF]. It suggests that NAO is related to a lack of cooling effect rather than high nasal resistance. Perhaps the truth lies in a combination of both [START_REF] Pawar | Objective measures in aesthetic and functional nasal surgery: perspectives on nasal form and function[END_REF].

Nasal valve area is a key-zone in the sensation of NAO [START_REF] Lindemann | Impact of menthol inhalation on nasal mucosal temperature and nasal patency[END_REF]. It is also difficult to treat surgically, hence the major interest of understanding its pathophysiology (48) [START_REF] Nyte | Hyaluronic acid spreader-graft injection for internal nasal valve collapse[END_REF]. Pre-and postoperative analyzes of patients with NVC or high and anterior SD show a link between a sudden total pressure drop and NAO perception (Figure 3a, b) (26) [START_REF] Kimbell | Changes in nasal airflow and heat transfer correlate with symptom improvement after surgery for nasal obstruction[END_REF]. In this area velocities are often the highest: Liu finds that the caudal septal deviations responsible for NVC are both correlated with higher velocities and NAO perception [START_REF] Sullivan | Perception of better nasal patency correlates with increased mucosal cooling after surgery for nasal obstruction[END_REF]. The study of pressures and velocities is easy and must be systematically undertaken when performing CFD.

Regarding WSS, it is in nasal valve area we find the most important friction areas [START_REF] Kimbell | Changes in nasal airflow and heat transfer correlate with symptom improvement after surgery for nasal obstruction[END_REF].

But the study of WSS also brings a lot in the analysis of patients with an ENS. Indeed, in this specific population, ENS6Q questionnaire was correlated with WSS: the absence of friction around the inferior turbinate causes the paradoxical NAO perception (50) [START_REF] Li | Computational fluid dynamics and trigeminal sensory examinations of empty nose syndrome patients[END_REF]. However, while presented results seem interesting, there is still a lack of data to interpret the WSS as well as possible. It will be necessary to correlate WSS and patient perception, before and after surgery.

Concerning the airway, several authors report a strong correlation between the patient's obstruction and the path taken by the air inside the nasal cavity. Patients with a decreased flow at the middle meatus appear more uncomfortable in case of SD [START_REF] Geißler | Rhinomanometry and acoustic rhinometry[END_REF].

Of the studies reviewed, two do not include validated assessment of NAO but simply report a clinical difference between the groups of patients [START_REF] Kimbell | Changes in nasal airflow and heat transfer correlate with symptom improvement after surgery for nasal obstruction[END_REF][START_REF] Zhao | Regional peak mucosal cooling predicts the perception of nasal patency[END_REF]. However, given the number of studies, we have incorporated all those that report the patient's perception. Moreover, because of the small number of patients included, sometimes correlations seem to be emerging but without satisfactory statistical test. It is therefore necessary to conduct other studies on a larger number of patients.

Other data remains to be integrated to make the simulations more and more precise: in the future, we have to be able to perform more often inspiratory and expiratory simulations, following a normal breathing (and not a constant incoming airflow), to allow a deformation of nasal walls which have a major importance regarding NVC, to integrate new data such as hygrometry or to vary the temperature of the nasal mucosa according to inspiration or expiration. This will allow more reliable extrapolation of CFD results.

CONCLUSION

The multiplication of CFD studies in the nose helps to better understand NAO. The clinical interpretation of previously unknown data, such as WSS or HF, opens up new horizons in the understanding of this symptom. Further studies on larger cohorts of patients need to be undertaken. However, we must be keep in mind that NAO may include a notion of subjectivity, non "CFDable".

  patients; -Pre vs post operative situations. All results are presented for the most obstructed side. SD = Septal Deviation. ITH = Inferior Turbinate Hypertrophy. NVC = Nasal Valve Collapse. Bold and * indicates a statistically significant correlation or difference (p<0.05). NS = Non significative test. If empty, not analyzed.

Figures

  Figures

Figure 1 .

 1 Figure 1. PRISMA Flow Diagram.

  Figure 1. PRISMA Flow Diagram.

Figure 2 .Figure 3 .

 23 Figure 2. Spatial distributions of Heat Fluxes on the nasal septum of a patient with nasal

Figure 3a .

 3a Figure 3a. shows the distributions of the relative total pressure during inspiration of the

Table 1 .

 1 Clinical assessment of NAO according to studies.

	Tables
	Article	Clinical assessment of NAO
		Subjective	Objective
	Casey et al. 2017	NOSE / VAS	/
	(19)
	Gaberino et al. 2017	NOSE / VAS	/
	(20)
	Hildebrandt et al.	Improvement after	/
	2013 (26)	surgery
	Kim et al. 2014 (22)	Patients with NAO vs	/
		healthy subjects
	Kimbell et al. 2012	NOSE / VAS	/
	(23)
	Kimbell et al. 2013	NOSE / VAS	/
	(21)
	Li, Farag, Leach et al.	NOSE / SNOT-22 / ENS6Q Rhinomanometry / Acoustic rhinometry
	2017 (30)
	Li, Farag, Maza et al.	NOSE / SNOT-22 / ENS6Q	/
	2017 (31)
	Liu et al. 2012 (29)	VAS	/
	Lu et al. 2014 (27)	/	Rhinomanometry / Acoustic rhinometry
	Rhee et al. 2012 (9)	NOSE / VAS	/
	Sullivan et al. 2013	NOSE / VAS	/
	(28)
	Zhao et al. 2013 (25)	VAS	Rhinomanometry

Table 2 .

 2 The table reports all CFD-calculated values. Comparison between CFD and:-NOSE and VAS scores; -Symptomatic and non symptomatic

Funding source: none Financial disclosure: none

Conflict of interest: none

Contrariwise, when the inferior turbinate is missing, there is an increase of the flow in the middle meatus [START_REF] Hsu | Evaluation of nasal patency by visual analogue scale/nasal obstruction symptom evaluation questionnaires and anterior active rhinomanometry after septoplasty: a retrospective one-year follow-up cohort study[END_REF][START_REF]Casey_Correlation between Subjective Nasal Patency and Intranasal Airflow Distribution[END_REF]. Changing the passage of air through the nasal cavity, even without deterioration of nasal resistance, changes the perception of patients. It is likely that this impairment is multifactorial: thanks to CFD, it appears that mucosal cooling effect, HF and humidity exchanges are related to the symptoms brought by the patients [START_REF] Otto | Computational fluid dynamics evaluation of posterior septectomy as a viable treatment option for large septal perforations[END_REF]. As such, we did not find study involving humidity exchange and NAO evaluation.

Li et al. and Lu et al. compare acoustic rhinometry with MCAs and found no significant

difference [START_REF] Hildebrandt | Evaluation of the intranasal flow field through computational fluid dynamics[END_REF][START_REF] Liu | Effects of septal deviation on the airflow characteristics: using computational fluid dynamics models[END_REF][START_REF] Li | Computational fluid dynamics and trigeminal sensory examinations of empty nose syndrome patients[END_REF]. This confirms that acoustic rhinotry is a reliable morphological examination. On the other hand, although often presented as functional, we think that it is an anatomical examination. It is therefore difficult to make correlations between subjective perceptions and cross section areas.

Limits

Performing comparisons between clinical data and CFD remains difficult. To do this, authors have several possibilities: they can either compare healthy to pathological subjects (for example a patient population with SD and NAO vs a population without SD and NAO), or compare the same patient before and after surgery. Then, they compare clinical tests with CFD. It remains difficult to carry out such studies on large populations of subjects, for many reasons. First of all, the realization of postoperative CT scan is not always justified, and it is necessary to obtain the approval of the local Ethics Committee, as CTscans exposes the patient to additional radiation. Moreover, the calculation and implementation times are still long and often require the presence of a computer engineer, programs available at the moment not being "user-frendly" for a surgeon. New software seems to be easier to use but we do not have the experience [START_REF] Balakin | Aerodynamic evaluation of the empty nose syndrome by means of computational fluid dynamics[END_REF].


Ulysse Faure

Csaba Domokos
email: csaba.domokos@milvus.ro

Agathe Leriche

Bogdan Cristescu

Brown bear den characteristics and selection in eastern Transylvania, Romania

Keywords: bear conservation, disturbance ecology, habitat ecology, MaxEnt, topography, Ursidae

Dens are important for species that need to survive and reproduce during harsh winters. Brown bears (Ursus arctos) in Romania, listed by the European Union as a population of concern, use dens for several months each year. To date, few quantitative assessments of denning habitat have been carried out for this population or others in Europe. In 2008-2013 and 2015-2017, we used local knowledge and telemetry data from brown bears fitted with GPS collars to identify 115 winter dens and eight open ground nests used by bears in eastern Transylvania, Romania. We located most dens in mountainous areas (64%) and fewer in foothills (36%). Den entrances in mountainous areas were significantly narrower than entrances in foothills, likely due to the need for reduced thermal loss during more severe winters at higher elevations. We selected seven habitat characteristics (abiotic and biotic) and human-related covariates associated with known locations of dens and open nests to identify potential brown bear denning habitat using maximum entropy modeling. We found that terrain ruggedness was the single most important factor when predicting bear denning habitat. The habitat map derived from this study can be used in the future to safeguard bear denning areas from potential human disturbances.

Some mammals respond to fluctuations in resource availability by migrating seasonally between wintering and breeding grounds [START_REF] Fryxell | Causes and consequences of migration by large herbivores[END_REF][START_REF] Alerstam | Long-distance migration: evolution and determinants[END_REF][START_REF] Avgar | On the adaptive benefits of mammal migration[END_REF], whereas others undergo prolonged hibernation or torpor (e.g., [START_REF] Millesi | Sex and age differences in mass, morphology, and annual cycle in European ground squirrels, Spermophilus citellus[END_REF][START_REF] Ortmann | Regulation of body temperature and energy requirements of hibernating alpine marmots (Marmota marmota)[END_REF][START_REF] Kryštufek | Glis glis (Rodentia: Gliridae)[END_REF]). Decreasing metabolism and sleeping inside a den for the winter allows some mammals to avoid the energetic demands of maintaining high body temperatures in winter, when food often is scarce [START_REF] Boyer | Molecular and metabolic aspects of mammalian hibernation[END_REF]. Winter dens also can provide a safe and sheltered environment for giving birth, nursing, and early cub growth, for some bear species (Ursidae- [START_REF] Oli | Denning ecology of black bears in the White River National Wildlife refuge, Arkansas[END_REF][START_REF] Seryodkin | Denning ecology of brown bears and Asiatic black bears in the Russian Far East[END_REF][START_REF] Manchi | Denning behaviour of Scandinavian brown bears Ursus arctos[END_REF][START_REF] Derocher | Sea ice and polar bear den ecology at Hopen Island, Svalbard[END_REF].

Brown bears (Ursus arctos) spend considerable time in their dens [START_REF] Schoen | Denning ecology of brown bears on Admiralty and Chichagof Islands[END_REF]), especially at northern latitudes [START_REF] Manchi | Denning behaviour of Scandinavian brown bears Ursus arctos[END_REF]. Winter denning is essential for female brown bears because they use the security of winter dens to give birth to blind cubs that have one of the lowest cub to mother mass ratio of all placental mammals: ~0.5 kg, or < 1% [START_REF] Huber | Food intake and mass gain of hand-reared brown bear cubs[END_REF][START_REF] Pasitschniak-Arts | Ursus arctos[END_REF]. During denning, bears do not eat, but instead rely on fat reserves accumulated in the predenning period [START_REF] Nelson | Behavior, biochemistry, and hibernation in black, grizzly, and polar bears[END_REF], resulting in the loss of 22-40% of their predenning body weight [START_REF] Kingsley | Patterns of weight gain and loss for grizzly bears in northern Canada[END_REF]. Bears also have adapted to winter denning by not drinking water, urinating, or defecating [START_REF] Linnell | How vulnerable are denning bears to disturbance?[END_REF], and reducing their metabolism by ~70% [START_REF] Watts | Energetic cost of winter dormancy in grizzly bear[END_REF], heart rate from 75 to 26 bpm [START_REF] Jørgensen | Low cardiac output as physiological phenomenon in hibernating, free-ranging Scandinavian brown bears (Ursus arctos) -an observational study[END_REF], and body temperature from 37°C to 33°C [START_REF] Hissa | Seasonal patterns in the physiology of the European brown bear (Ursus arctos arctos) in Finland[END_REF].

Flooding, especially following the melting of snow in spring [START_REF] Schoen | Denning ecology of brown bears on Admiralty and Chichagof Islands[END_REF][START_REF] Huber | Denning of brown bears in Croatia[END_REF] and den collapse [START_REF] Miller | Denning ecology of brown bears in southcentral Alaska and comparisons with a sympatric black bear population[END_REF]) sometimes can affect denning bears. Human disturbances, such as hunting and resource extraction, also constitute threats that in some cases result in bears abandoning dens [START_REF] Swenson | Winter den abandonment by brown bears Ursus arctos: causes and consequences[END_REF]. Den abandonment can have potentially lethal effects on litters of cubs [START_REF] Linnell | How vulnerable are denning bears to disturbance?[END_REF]. Consequently, the selection of a denning site with a low probability of human disturbance is important for winter survival, successful reproduction, and safety from humans [START_REF] Sahlén | Den entry behavior in Scandinavian brown bears: implications for preventing human injuries[END_REF]. From a management standpoint, educating the public to minimize human-bear interactions during winter denning may promote social acceptance of brown bears, thereby assisting conservation efforts [START_REF] Sahlén | Den entry behavior in Scandinavian brown bears: implications for preventing human injuries[END_REF].

Previous studies have documented several den types used by brown bears, including natural rock cavities, excavations in or under live trees or snags, excavations in soil, and surface beds [START_REF] Schoen | Denning ecology of brown bears on Admiralty and Chichagof Islands[END_REF][START_REF] Linnell | How vulnerable are denning bears to disturbance?[END_REF][START_REF] Seryodkin | Denning ecology of brown bears and Asiatic black bears in the Russian Far East[END_REF]. Dens typically consist of a porch, entrance, tunnel, and bedding chamber [START_REF] Miller | Denning ecology of brown bears in southcentral Alaska and comparisons with a sympatric black bear population[END_REF]; Fig. 1). The porch is defined as the front portion of the den and often is constructed from soil excavated from inside the den. Bears occasionally build a bed on the porch for sunbathing during winter [START_REF] Harding | Den-site characteristics of arctic coastal grizzly bears (Ursus arctos L.) on Richards Island, Northwest Territories, Canada[END_REF][START_REF] Étienne | L'Ours brun d'Europe, biologie et histoire, des Pyrénées à l'Oural. Biotope, Mèze (Collection Parthénope)[END_REF]. Within the bedding chamber, bears typically build beds using branches, leaves, grass, and other materials, available in the immediate vicinity [START_REF] Schoen | Denning ecology of brown bears on Admiralty and Chichagof Islands[END_REF][START_REF] Li | Characteristics of dens and selection of denning habitat for bears in the South Xiaoxinganling Mountains, China[END_REF].

Past studies have identified a variety of factors responsible for brown bears selecting den sites, including land cover, food resources, topography, and anthropogenic features [START_REF] Petram | Human influence on the choice of winter dens by European brown bears in Slovenia[END_REF][START_REF] Elfström | Selection of denning habitats by Scandinavian brown bears Ursus arctos[END_REF][START_REF] Libal | Microscale den-site selection of grizzly bears in southwestern Yukon[END_REF][START_REF] Pigeon | Den selection by grizzly bears on a managed landscape[END_REF]. The den site selection process occurs at multiple scales. For example, at a broad scale, bears may select areas for denning that are located in alpine or old growth forests [START_REF] Schoen | Denning ecology of brown bears on Admiralty and Chichagof Islands[END_REF], and that have abundant and high-quality spring foods in the vicinity [START_REF] Pigeon | Den selection by grizzly bears on a managed landscape[END_REF]. At a more local scale, the presence of stabilizing den material such as large rocks [START_REF] Ciarniello | Denning behavior and den site selection of grizzly bears along the Parsnip River, British Columbia, Canada[END_REF] and horizontal cover at 1 m (approximate height of walking bears- [START_REF] Libal | Microscale den-site selection of grizzly bears in southwestern Yukon[END_REF]) positively influence den site selection.

Terrain also is a major factor that contributes to the process of den site selection [START_REF] Petram | Human influence on the choice of winter dens by European brown bears in Slovenia[END_REF][START_REF] Pigeon | Den selection by grizzly bears on a managed landscape[END_REF], with brown bears in some regions selecting only rugged terrain for denning [START_REF] Whiteman | Identifying key denning habitat to conserve brown bear (Ursus arctos) in Croatia[END_REF]. Past studies found that extreme ruggedness provides brown bears safety from human disturbance, because such terrain is difficult to access during winter months [START_REF] Nielsen | Modelling the spatial distribution of human-caused grizzly bear mortalities in the Central Rockies ecosystem of Canada[END_REF]. Forest cover also provides bears food and refuge to evade people [START_REF] Pop | Combining resource selection functions and home-range data to identify habitat conservation priorities for brown bears[END_REF]. In some regions, brown bears select certain slopes to maximize the buildup of snowfall and thermal insulation in the den [START_REF] Craighead | Data on grizzly bear denning activities and behavior obtained by using wildlife telemetry[END_REF][START_REF] Vroom | The ecology of winter den sites of grizzly bears in Banff National Park, Alberta[END_REF][START_REF] Schoen | Denning ecology of brown bears on Admiralty and Chichagof Islands[END_REF]. Distance to roads and agricultural land can have a negative impact on site selection, possibly because brown bear associate roads, infrastructure, and other developed land, with people [START_REF] Swenson | Winter den abandonment by brown bears Ursus arctos: causes and consequences[END_REF][START_REF] Sahlén | Brown bear den site concealment in relation to human activity in Sweden[END_REF].

To date, most research on brown bear den site selection has been carried out in North America [START_REF] Ciarniello | Grizzly bear habitat selection is scale dependent[END_REF][START_REF] Goldstein | Brown bear den habitat and winter recreation in south-central Alaska[END_REF][START_REF] Hodder | Application of a species distribution model to identify and manage bear den habitat in central British Columbia, Canada[END_REF][START_REF] Pigeon | Den selection by grizzly bears on a managed landscape[END_REF][START_REF] Pigeon | Assessing den selection and den characteristics of grizzly bears[END_REF]; few studies in Europe have described brown bear dens and the process of den site selection [START_REF] Huber | Denning of brown bears in Croatia[END_REF][START_REF] Petram | Human influence on the choice of winter dens by European brown bears in Slovenia[END_REF][START_REF] Elfström | Selection of denning habitats by Scandinavian brown bears Ursus arctos[END_REF]. Despite Romania hosting the largest brown bear population in Europe (> 6,000 individuals), excluding European Russia [START_REF] Swenson | Action plan for the conservation of the brown bear in Europe[END_REF][START_REF] Chapron | Recovery of large carnivores in Europe's modern human-dominated landscapes[END_REF], denning of brown bear in Romania has received little attention from researchers over the past 50 years [START_REF] Almășan | Zur Kenntnis des Rumänischen Karpatenbär[END_REF].

The purpose of this study was to describe brown bear dens and gain an understanding of the factors responsible for selection of den sites by bears in eastern Transylvania, Romania. First, we compared the structure and dimensions of known dens found in the foothills and mountains. We then used these data to test the structural thermal insulation hypothesis, which states that bears construct dens with longer tunnels and smaller entrances in mountainous regions presumably because winters are more severe at higher elevations. We also tested several den site selection hypotheses using data from known den locations: i) the risk minimization hypothesis, which states that terrain ruggedness and forest cover are positively associated with den sites; ii) the snowpack thermal insulation hypothesis, which claims that dens will occur where snow cover is more contiguous, deeper, and longer lasting; and iii) the human avoidance hypothesis, which states that high road density, infrastructure, and human-modified habitats are negatively associated with den sites. Lastly, we used habitat data for confirmed dens and Species Distribution Models (SDMs) to predict the location of sites suitable for denning within the study area.

Materials and Methods

Study area.-The study was conducted in eastern Transylvania, Romania (Fig. 2). The east section of the study area (1,121 km 2 ; hereafter, mountains) is mountainous and includes part of Romania's eastern Carpathians (elevation range: 515-1,783 m). More than a third of the eastern study area (35.6%) is protected by the European Union's (EU) Habitats Directive as a Natura 2000 site. Forests account for 68% of total land cover and include: European beech (Fagus sylvatica; 25%); Norway spruce (Picea abies), silver fir (Abies alba), and European larch (Larix decidua; jointly 26%); and mixed forests (17%). Approximately 20% of the land is pastureland used to raise livestock from late spring to early autumn. A small portion of the study area (0.5%) also includes human settlements, weekend homes, an active stone quarry, and a few small ski resorts (CORINE Land Cover database 2012). The mountainous area has a road density (mean length of roads used by motorized vehicles per km 2 ) of 0.62 km/km 2 . The vast majority of the transportation infrastructure consists of gated roads maintained by local forestry services. Other than a few paved national and county roads and roads leading to forest parcels with active logging operations, the roadways are not cleared of snow during the denning period, making them inaccessible to most vehicles.

The west section of the study area consists mainly of the foothills of the eastern Carpathian Mountains (1,999 km 2 ; hereafter, foothills), extending west of the Gurghiu Mountains and south of the Gurghiu River (elevation range: 310-1,090 m). A small proportion of the area (< 20%) is part of the Natura 2000 European network of protected areas. The foothills are largely accessible to humans, hosting numerous settlements (7% of the study area-CORINE Land Cover database 2012), ranging from small villages to towns. The main land uses (60%) are arable and pastoral farming. Corn, alfalfa, oat, potato, and wheat croplands occur mainly at lower elevations, and seminatural pastures, hayfields, and orchards at higher elevations. Forest succession occurs in abandoned fields where pioneer species, including blackthorn (Prunus spinosa), common hawthorn (Crataegus monogyna), and common sea buckthorn (Hippophae rhamnoides), form almost impenetrable thickets that develop into deciduous forests (27% of the area) of European hornbeam (Carpinus betulus), oak (Quercus sp.), European beech, and black locust (Robinia pseudoacacia). Plantations of Scots pine (Pinus sylvestris; 0.3%) and mixed forests (0.2%) also are present. Forests cover much of the high elevations and narrow valleys, both of which generally are difficult to access by people. The foothills have a road density of 1.2 km/km 2 . Major roads are cleared of snow throughout the denning period. Most forestry roads are not gated and, because winters are relatively mild, some roads are periodically accessible to vehicles throughout winter.

Identifying dens.-We located bear dens in 2008-2013 and 2015-2017. In the mountains, we found dens with help from foresters and game managers. In the foothills, we located dens with assistance from game managers, pastoralists, forest owners, mushroom hunters, and berry pickers. In 2015-2017, we identified some dens based on telemetry data from bears fitted with GPS collars. We captured bears in cage traps and immobilized them using a dosage of 35 µg Medetomidine and 2-5 mg Zoletil (tiletamine-zolazepam) per estimated body mass (kg) using dart syringes fired from a Dan-Inject CO 2 Injection Rifle (Model J.M.SP; Dan-Inject ApS, Kolding, Denmark). We fitted captured bears with Vectronic Aerospace GPS Pro Light-4 GPS-GSM collars (Vectronic Aerospace, Berlin, Germany). We set collars to record GPS locations of bears every 1 h, including the denning period.

Unless we were absolutely certain that the structure was unoccupied, den visits only occurred outside the denning period, so as to avoid disturbance and prevent potential den abandonment. We carefully inspected each suspected denning structure and identified it as a den only if it contained clearly visible claw marks (e.g., from the digging process, or from usage), bear hair, or clearly recognizable beds. We recorded every den's location, elevation (as determined with a Garmin GPSMAP 60CSx GPS unit; Garmin International, Olathe, Kansas), and land cover type in the immediate vicinity (10-m radius). We only recorded the location of aboveground open nests (excavated, circular, shallow depressions usually containing bedding material) used by wintering bears.

Den characterization.-We described the dimensions of dens (excluding open nests) using the following measurements: height and width of entrance; height, width, and length of tunnel and chamber; and porch length. All measurements were collected by the same person using a measuring tape (cm).

The structural thermal insulation hypothesis predicted that bears would construct dens with longer tunnels and smaller entrances in mountainous regions because winters are more severe. We therefore compared the measurements of dens between the foothills and mountains. Because the data were not normally distributed according to Shapiro-Wilk tests and graphical investigations, we carried out nonparametric Wilcoxon rank sum tests to contrast the den structural metrics. We used R core functions and the dplyr and ggpubr packages in R studio (v 1.1.447-R Development Core Team 2017), and assessed the influence of outliers on our results.

Modeling of Denning Habitat

Modeling method.-Species Distribution Modeling can be used to help predict the distribution or habitat of a species, such as habitat used for denning as we do here, across geographic space and time using environmental data. In brief, we used maximum entropy modeling (MaxEnt) to assess the probability of presence of brown bear dens across the study area based on habitat data collected at known den locations. MaxEnt identifies denning habitat by finding the distribution of maximum entropy (i.e., closest to uniform- [START_REF] Phillips | Maximum entropy modeling of species geographic distributions[END_REF][START_REF] Elith | A statistical explanation of MaxEnt for ecologists[END_REF]; where environmental covariates in the predicted extent matched average values recorded at confirmed denning locations. MaxEnt has been shown to be well-adapted to presence-only data and has consistently demonstrated performance that is competitive with other methods, such as generalized linear models, generalized additive models, Genetic Algorithms for Rule-set Prediction (GARP), or bioclimatic envelopes (BIOCLIM- [START_REF] Elith | Novel methods improve prediction of species' distributions from occurrence data[END_REF][START_REF] Phillips | Maximum entropy modeling of species geographic distributions[END_REF].

We fitted models by investigating linear, quadratic, and product feature types and all combinations of those features. The machine learning fitting process was set to end when training gain fell below the threshold value of 0.0001, or after 1,000 iterations. The background data set, which is a random sample of points from the study area [START_REF] Elith | A statistical explanation of MaxEnt for ecologists[END_REF], was set with 10,000 maximum samples within a restricted area defined as the Minimum Convex Polygon of the den occurrences [START_REF] Webber | Modelling horses for novel climate courses: insights from projecting potential distributions of native and alien Australian acacias with correlative and mechanistic models[END_REF]). The regularization multiplier was set to 1 to minimize the problem of model overfitting, while not producing a widely generalized model [START_REF] Phillips | Maximum entropy modeling of species geographic distributions[END_REF]. A jackknife test also was performed to obtain an estimate of individual variable importance [START_REF] Elith | A statistical explanation of MaxEnt for ecologists[END_REF].

Biased samples can be problematic in distribution modeling. The data collection methods for the den occurrence records could not assure a nonbiased data set. Most dens were discovered opportunistically by people accessing areas throughout the region, including wild boar hunters, shepherds, foresters, mushroom hunters, and berry pickers. At least some of these people likely accessed areas from a nearest road, which could have resulted in spatial distribution of dens concentrated around roadways. To reduce potential spatial autocorrelation caused by proximity to roads, we introduced a bias correction to the models, using background manipulation as described in past studies [START_REF] Phillips | Sample selection bias and presence-only distribution models: implications for background and pseudoabsence data[END_REF][START_REF] Kramer-Schadt | The importance of correcting for sampling bias in MaxEnt species distribution models[END_REF]). We calculated the Euclidean distance to forestry roads, then reversed the distance values and included the final raster as a bias file in MaxEnt.

We fitted six models (features L, Q, LP, LQ, QP, LQP) with and without bias correction, totaling 12 models. We then produced 10 replicates for each feature and condition combination (k-fold cross-validation method). For modeling purposes, we reduced the overall data set of 123 known dens to 119 dens. This ensured that each modeling unit (300 m × 300 m) contained one den location only.

Explanatory covariates.-We extracted spatial covariates using ArcGIS (v10.3.1-ESRI 2015) and QGIS software (v2.0.1-QGIS Development Team 2013) at a resolution of a modeling unit (300 m × 300 m). We included as covariates abiotic and biotic habitat characteristics, as well as human-related features. The risk minimization hypothesis predicts that bears should den in rugged terrain and forested areas; we therefore estimated terrain ruggedness (accessibility/disturbance) and forest cover (foraging habitat and refuge). The snowpack thermal insulation hypothesis predicts that dens occur where snow cover is most extreme; to index snow occurrence, we calculated aspect (thermal insulation), assuming that northern and western aspects retain most snow. The human avoidance hypothesis predicts that denning takes place in areas with low densities of roads and infrastructure, and where disturbed habitats are in low occurrence; accordingly, we generated density of roads, percentage of artificial area, and percentage of humangenerated open fields (disturbances).

We estimated ruggedness (TRugg) using a terrain ruggedness index (TRI) calculated from a 30-m Digital Elevation Model (DEM) from the GMES RDA project (EU-DEM, https://www. eea.europa.eu/; Supplementary Data SD1). The TRI provides a quantitative measure of topographic heterogeneity, calculating the sum change in elevation between a grid cell and its eight neighboring cells [START_REF] Riley | A terrain ruggedness index that quantifies topographic heterogeneity[END_REF].

The CORINE Land Cover 2012 polygonal vector layer (CORINE Land Cover database 2012) was used to delineate forest cover percentage (ForestC; Supplementary Data SD2). Due to sample size limitations, forest cover pooled conifer, deciduous, and mixed, forests. We undertook calculations with Focal Statistics within a 495-m rectangle buffer around each cell [START_REF] Elfström | Selection of denning habitats by Scandinavian brown bears Ursus arctos[END_REF].

We derived den aspect from the GMES RDA DEM. Because aspect is a circular variable with values ranging from 1 to 360, it was converted to two continuous aspect-index variables (Northness, NorthI [Supplementary Data SD3] and Eastness, EastI [Supplementary Data SD4]), using the "raster" and "sp" packages in RStudio v0.99.4 (R Development Core Team 2017). These aspect indices range from -1 (southern/western aspects) to 1 (northern/eastern aspects).

We exported the linear vector layer of the Romanian road network from the OpenStreetMap (OSM) project database (http://www.openstreetmap.org/; Open Database License [ODbl] v1.0) and used it to calculate the density of paved roads (per 100 km 2 ) within a 500-m circular buffer [START_REF] Elfström | Selection of denning habitats by Scandinavian brown bears Ursus arctos[END_REF]. The covariate (RoadD) included all paved OSM road types, which are used in the study area by cars, horse carts, and pedestrians (Supplementary Data SD5).

We obtained the percentage of artificial area (human infrastructure, ArtifC) and percentage of human-generated open fields (crops and pastures, OpenC) from the CORINE Land Cover 2012 polygonal vector layer (CORINE Land Cover database 2012). We used Focal Statistics within a 495-m rectangle buffer around each cell to calculate these percentages [START_REF] Elfström | Selection of denning habitats by Scandinavian brown bears Ursus arctos[END_REF].

We inspected correlations between continuous covariates using ENMTools software (ENMTools 1.4.4 version- [START_REF] Phillips | Maximum entropy modeling of species geographic distributions[END_REF][START_REF] Warren | Environmental niche equivalency versus conservatism: quantitative approaches to niche evolution[END_REF]) and Pearson's correlation tests with 95% confidence intervals. We did not include variables in the same model if correlated r > |0.65| [START_REF] Kuemmerle | European bison habitat in the Carpathian Mountains[END_REF]). We found ForestC and OpenC were negatively correlated (r = -0.82; Supplementary Data SD6) and RoadD was positively correlated with ArtifC (r = 0.71). For easy interpretation and straightforward applicability to management, ForestC and RoadD were used to model the distribution of potential denning areas and OpenC and ArtifC were excluded.

Model selection and projection.-We used the Akaike Information Criterion corrected for small sample sizes (AICc) to select the top model of 12. We calculated the AICc value for each replicate of each feature combination using ENMTools software [START_REF] Hurvich | Model selection for extended quasi-likelihood models in small samples[END_REF][START_REF] Warren | Ecological niche modeling in Maxent: the importance of model complexity and the performance of model selection criteria[END_REF]. The most parsimonious model [START_REF] Mazerolle | Making sense out of Akaike's Information Criterion (AIC): its use and interpretation in model selection and inference from ecological data[END_REF]) was projected over the full study area, with the same parametrization as the one used during the fitting process, excepting the output format that was set to logistic with 10 replicates. The Area Under the receiver operating Curve function (AUC) was used to assess the performance of our top model. The AUC represents the probability of a random presence location to be ranked with a higher value than a random background sample [START_REF] Fielding | A review of methods for the assessment of prediction errors in conservation presence/absence models[END_REF][START_REF] Merow | A practical guide to MaxEnt for modeling species' distributions: what it does, and why inputs and settings matter[END_REF]. We also investigated the degree of model extrapolation by generating a Multivariate Environmental Similarity Surface (MESS; [START_REF] Elith | The art of modelling range-shifting species[END_REF]) for brown bear denning, which was spatially constrained to the study area extent. MESS analysis quantified the similarity in covariates between any given location in the projection data set and the locations in the reference (training) data set [START_REF] Elith | The art of modelling range-shifting species[END_REF].

Results

Den locations.-We identified a total of 115 bear dens and eight open nests. We located 74 dens (60%) and three open nests (2%) in the mountains, and 41 dens (33%) and five open nests (4%) in the foothills. Overall, 10 dens and one open nest were located using telemetry. With the exception of one open nest in the mountains used by a collared male, all open nests were used by females with cubs (n = 7), which was confirmed visually by the authors or locals.

Den characterization.-We found that dens were completely dry when inspected, with the exception of four dens in the foothills that were flooded. Most mountain dens (73%; n = 54) were excavated (31 under the roots of trees; 23 under rocks), whereas 27% (n = 20) were natural cavities located in or under rocks or between large boulders. We identified beds in 80% of the dens (n = 59); many dens contained bedding material (n = 37). We also documented tunnels in 54% (n = 40) of the structures. Elevations of mountain den locations ranged from 486 to 1,607 m (x = 1,118 m; SD = 249 m) in different types of forests: 32% (n = 24) in deciduous forests, 38% (n = 28) in mixed forests, and 30% (n = 22) in conifer forests.

A larger proportion of dens in the foothills than in the mountains (95%; n = 39) were excavated (24 were dug under the roots of trees, and 15 were dug into the sides of slopes); the remaining two dens in the foothills included one natural cavity in a rock, and one under the roots of a tree. Tunnels were present in 27% of the structures (n = 11). We identified beds in 44% of the dens (n = 18), and some of the dens (n = 9) contained bedding material. Elevations of foothills den locations ranged from 393 to 884 m (x = 491 m; SD = 77 m) in deciduous forests (68%; n = 28), thickets on abandoned agricultural land (29%; n = 12), and a planted patch of conifer forest (2%; n = 1). Elevations of some foothills dens ranged higher than the low range of elevations of mountain dens, because foothills and mountains were delineated by biogeographical region rather than exclusively elevation.

We found some support for the structural thermal insulation hypothesis. Of the nine den structural metrics (Table 2), only entrance width differed significantly between mountains and foothills. Entrances of mountain dens were significantly narrower than those in the foothills (P = 0.004). We only discussed our results derived from our full data set because results were similar with and without outlying values (Supplementary Data SD7).

Habitat modeling of den sites.-The top two models in the candidate set included linear and quadratic features (LQ) without and with a bias file, respectively (mean AUC = 0.77 ± 0.05; mean AUC = 0.78 ± 0.06; Supplementary Data SD8). The LQ model without the bias file had the lowest AICc score (Supplementary Data SD8) and low standard deviation (Supplementary Data SD9). No other models in the candidate set had an ∆AICc < 2.

Our prediction that bears would den in rugged terrain and forested areas was partially supported. Terrain ruggedness (TRugg) was the best predictor of den presence (93.0% of model contribution; Fig. 3). The top model predicted a 0.5 probability of den occurrences at a TRugg index value of 21.7, and a maximum probability of presence at a TRugg of 41.6 (Fig. 3). The slight decrease in probability of presence with very high TRugg suggests that very steep areas are less suitable den sites (Fig. 3). In contrast, forest cover had negligible effect on model predictions (0.8%) (Supplementary Data SD10).

We also found partial support for our prediction that denning would occur where snow cover is most extreme. Aspect contributed relatively little to model predictions (5.9% in total) with slightly higher probabilities on northern-and westernfaced slopes (Supplementary Data SD10).

The prediction that bears would den at low densities of human disturbance as indexed by main transportation infrastructure received little support. Only 0.2% of model predictions were explained by paved-road density (RoadD; Supplementary Data SD10).

Multivariate Environmental Similarity Surface (MESS) results indicate that there were no overextrapolations of the top model when predicting the presence of dens across the entire area (Supplementary Data SD11). Our model predicted that rugged areas in the Northeast and Southwest sections of the study area were most suitable for brown bear denning and lowlands, plateaus, and river floodplains were the least suitable (Fig. 4).

Discussion

Habitat modeling aids in landscape-level planning, which in turn can guide the development of conservation and management strategies. Efforts to model brown bear habitat selection in Eastern Europe have lagged compared to similar work ongoing elsewhere in Europe and North America. In eastern Transylvania, Romania, we found that bears mostly excavated their dens (76%) in rugged terrain, which generally has low levels of human activity in winter compared to flatter regions. Areas with high potential for extreme snow conditions, or human disturbance associated with paved roads, may have little influence on the process of den site selection under current study area conditions.

All types of dens (excavations, caves, and nests) that we documented in this study have been described in previous studies (e.g., [START_REF] Ciarniello | Denning behavior and den site selection of grizzly bears along the Parsnip River, British Columbia, Canada[END_REF]. Studies in the United States (most of Alaska; Kodiak Island; Montana; Rocky Mountains and Yellowstone National Park), Canada (Banff and Jasper National Parks; the Yukon), and Sweden [START_REF] Linnell | How vulnerable are denning bears to disturbance?[END_REF] also found that most dens used by brown bears were excavated. We found that excavated dens were dug into the sides of slopes and often under the roots of trees or under rocks, which presumably decreases the risk of den collapse [START_REF] Ciarniello | Denning behavior and den site selection of grizzly bears along the Parsnip River, British Columbia, Canada[END_REF]. Similar to past studies in the United States (Yellowstone National Park, Admiralty Island, and Chichagof Island- [START_REF] Judd | Denning of grizzly bears in the Yellowstone National Park area[END_REF][START_REF] Schoen | Denning ecology of brown bears on Admiralty and Chichagof Islands[END_REF], Canada [START_REF] Vroom | The ecology of winter den sites of grizzly bears in Banff National Park, Alberta[END_REF][START_REF] Ciarniello | Denning behavior and den site selection of grizzly bears along the Parsnip River, British Columbia, Canada[END_REF], Croatia [START_REF] Huber | Denning of brown bears in Croatia[END_REF], and Italy [START_REF] Groff | Selection of denning caves by brown bears in Trentino, Italy[END_REF], we also found that bears used caves in rocks to den for the winter. Denning in caves was much more common in the mountains (20 of 22 caves) than in the foothills (2 of 22). All but one of the natural cavities used by bears were located in or under rocks, or between large boulders. The high frequency of rock formations in the volcanic Gurghiu Mountains [START_REF] Mason | Low-pressure evolution of magmas from the Calimani, Gurghiu and Harghita Mountains, East Carpathians[END_REF] likely explains the frequent presence of dens excavated under rocks and use of natural cavities in the mountains. A scarcity of rock formations in the foothills also could explain the lack of these den structures in this lower elevation region. Unlike excavated dens, natural cavities do not require energetic investment but may be less insulated.

We located only eight aboveground open nests. In Sweden and Norway, open nests primarily were used by males [START_REF] Manchi | Denning behaviour of Scandinavian brown bears Ursus arctos[END_REF][START_REF] Elfström | Effects of sex and age on den site use by Scandinavian brown bears[END_REF]. [START_REF] Manchi | Denning behaviour of Scandinavian brown bears Ursus arctos[END_REF] found that males that used open nests spent shorter time sleeping than males in other types of dens. The authors concluded that males might be able to use open nests both because of greater fat storage and less weight loss during winter denning compared to females. The need to shelter offspring from environmental conditions also might explain why females did not use open nests in Scandinavia. These open structures probably are less efficient than other den types in terms of reducing energy loss, particularly during periods when there is scant snow cover [START_REF] Elfström | Effects of sex and age on den site use by Scandinavian brown bears[END_REF]. In our study, seven of the eight open nests were used by females that had just given birth, and only one open nest was used by a male. One possible reason for the high use of open nests by females with cubs of the year in our study compared to [START_REF] Manchi | Denning behaviour of Scandinavian brown bears Ursus arctos[END_REF] was that winters in Romania are less severe than in Sweden.

Bedding material also contributes to thermal insulation in dens [START_REF] Craighead | Data on grizzly bear denning activities and behavior obtained by using wildlife telemetry[END_REF][START_REF] Mcloughlin | Denning ecology of barren-ground grizzly bears in the central Arctic[END_REF]. Although it is no surprise that we encountered more dens with bedding material in the mountains than in the foothills (50% of the dens versus 22%), we found fewer dens with bedding materials when compared to Croatia (86%-Huber and Roth 1997) and the United States (73%-Servheen and Klaver 1983; 76%- [START_REF] Judd | Denning of grizzly bears in the Yellowstone National Park area[END_REF]76%-Miller 1990). The site level availability of bedding materials might have influenced these findings, but was not quantified in our study.

While den entrance heights and tunnels did not vary between mountains and foothills, entrance widths differed, providing partial support for the structural thermal insulation hypothesis. Mountain den entrances (72.9 ± 29.5 cm) were approximately 10 cm narrower than those in the foothills (82.5 ± 20.3 cm). Results from past studies corroborate these findings, suggesting that bears do not adjust the height of the den entrance, as it is directly proportional to the size of the denning bear [START_REF] Servheen | Grizzly bear dens and denning activity in the Mission and Rattlesnake Mountains, Montana[END_REF], but instead they modify the width. Given similar den entrance heights, a narrower entrance is more easily obstructed by snow, effectively sealing off the den from the outside weather conditions and increasing thermal insulation (e.g., [START_REF] Craighead | Data on grizzly bear denning activities and behavior obtained by using wildlife telemetry[END_REF]. Even without a protective snow cover, a narrow entrance likely reduces heat loss [START_REF] Zhang | Factors predicting den use by maternal giant pandas[END_REF]).

We also found partial evidence toward the risk minimization hypothesis. Similar to other studies (e.g., [START_REF] Sahlén | Brown bear den site concealment in relation to human activity in Sweden[END_REF][START_REF] Smereka | Den selection by barren-ground grizzly bears, Mackenzie Delta, Northwest Territories[END_REF][START_REF] Whiteman | Identifying key denning habitat to conserve brown bear (Ursus arctos) in Croatia[END_REF], we identified terrain ruggedness as the most important variable for predicting locations where bears would den for the winter. When in proximity of human populated areas, bears may minimize the risk of encounters with humans by using areas of difficult access, including rugged terrain [START_REF] Martin | Coping with human disturbance: spatial and temporal tactics of the brown bear (Ursus arctos)[END_REF]. Selection of rugged terrain for den sites has also been reported for wolverines (Gulo gulo- [START_REF] May | Habitat characteristics associated with wolverine den sites in Norwegian multiple-use landscapes[END_REF] was found to be related to the presence of karstic formations used by denning bears in Croatia [START_REF] Whiteman | Identifying key denning habitat to conserve brown bear (Ursus arctos) in Croatia[END_REF]. In our study area, brown bears likely prefer rugged terrain because it offers greater availability of potential den sites and potentially because of lower levels of human disturbance in rugged areas [START_REF] Nielsen | Modelling the spatial distribution of human-caused grizzly bear mortalities in the Central Rockies ecosystem of Canada[END_REF][START_REF] May | Spatio-temporal ranging behaviour and its relevance to foraging strategies in wide-ranging wolverines[END_REF]. Terrain ruggedness at bear den sites was demonstrated to increase with decreasing distance to human presence and activity, such as roads and settlements [START_REF] Sahlén | Brown bear den site concealment in relation to human activity in Sweden[END_REF].

We found unexpected results regarding the relationship between forest cover and suitability for denning. Although all mountain dens and most foothill dens (70%) were located in forests, our model does not support preferential selection of forest habitat for denning (0.8%). Bears that overwinter close to human settlements and activity (as is the case in the foothills region of the study area) need to be concealed, reducing their detectability at den sites [START_REF] Sahlén | Brown bear den site concealment in relation to human activity in Sweden[END_REF]. Correspondingly, forest cover provides concealment, which is why forests may be used for denning by many large mammal species, such as Canada lynx (Lynx canadensis- [START_REF] Squires | Hierarchical den selection of Canada lynx in western Montana[END_REF]) and gray wolves [START_REF] Norris | Forest composition around wolf (Canis lupus) dens in eastern Algonquin Provincial Park, Ontario[END_REF][START_REF] Theuerkauf | Selection of den, rendezvous, and resting sites by wolves in the Bialowieza Forest, Poland[END_REF]. The surprising lack of association between forest cover and denning in the model results could be due to correlation between forest cover and ruggedness (r = 0.33). Both in the mountains and foothills, rugged areas are predominantly covered with forests. In contrast, woodless, open spaces are located in mostly flat terrain that can be grazed by sheep or cattle, and may be used for haymaking or as arable land. The disproportionately heavy weight of the ruggedness covariate in the model would mask the influence of forest cover, particularly when using a machine-learning algorithm.

We documented partial support for the snowpack thermal insulation hypothesis. The probability of den presence was highest on north-and west-facing slopes. However, the aspect covariate explained only a small percentage of the variability in den distribution in our study area. Previous studies have reported den site associations with all aspects: West [START_REF] Vroom | The ecology of winter den sites of grizzly bears in Banff National Park, Alberta[END_REF]; North [START_REF] Craighead | Data on grizzly bear denning activities and behavior obtained by using wildlife telemetry[END_REF][START_REF] Judd | Denning of grizzly bears in the Yellowstone National Park area[END_REF]); East [START_REF] Elfström | Selection of denning habitats by Scandinavian brown bears Ursus arctos[END_REF]; and South [START_REF] Schoen | Denning ecology of brown bears on Admiralty and Chichagof Islands[END_REF][START_REF] Miller | Denning ecology of brown bears in southcentral Alaska and comparisons with a sympatric black bear population[END_REF][START_REF] Groff | Selection of denning caves by brown bears in Trentino, Italy[END_REF]. Other authors found no preference for a specific aspect at den sites [START_REF] Huber | Denning of brown bears in Croatia[END_REF][START_REF] Seryodkin | Denning ecology of brown bears and Asiatic black bears in the Russian Far East[END_REF]. A review study concluded that aspect selection likely is influenced by the local stability of snow conditions [START_REF] Linnell | How vulnerable are denning bears to disturbance?[END_REF]. Several studies suggested that bears prefer slope exposures that accumulate the thickest snow cover, insulating the den [START_REF] Craighead | Data on grizzly bear denning activities and behavior obtained by using wildlife telemetry[END_REF]; Schoen et al. 1987). In our study area, west-facing slopes accumulate deeper snow cover due to predominant westerly winds that bring high precipitation, which discharges when encountering the west slopes; whereas slopes with northerly exposures retain snow for longer. The human avoidance hypothesis as indexed by main transportation infrastructure received scant support. Paved roads did not appear to have an effect on den site selection in our study area. [START_REF] Linnell | How vulnerable are denning bears to disturbance?[END_REF] suggested that bears den > 1 km from roads, whereas [START_REF] Elfström | Effects of sex and age on den site use by Scandinavian brown bears[END_REF] found that abandoned dens were closer to roads cleared of snow in the winter than successfully used dens. The authors concluded that plowed roads act as a source of disturbance, permitting increased human access to denning areas. [START_REF] Elfström | Selection of denning habitats by Scandinavian brown bears Ursus arctos[END_REF] documented that denning Scandinavian brown bears avoided roads with high disturbance potential, but denned near small roads that were not cleared of snow during winter, and large highways without parking possibilities. We hypothesize that the lack of effect of roads on denning habitat suitability might be due to most roads in our study area being either inaccessible or unused during the denning period. Similarly, [START_REF] Petram | Human influence on the choice of winter dens by European brown bears in Slovenia[END_REF] and [START_REF] Whiteman | Identifying key denning habitat to conserve brown bear (Ursus arctos) in Croatia[END_REF] found that distance to the nearest road had little effect on bear den site selection in Slovenia and Croatia.

Given the complex array of possible human disturbances to bear denning, the map of denning habitat suitability generated in this study is a useful conservation tool to assist in the spatial and temporal planning of human activities in areas exhibiting high probability of den presence, particularly during the denning period. We recommend that forestry management plans include provisions to minimize the amount and duration of forestry operations, or to avoid interventions altogether (road closures- [START_REF] Pigeon | Assessing den selection and den characteristics of grizzly bears[END_REF] in areas highly suitable for denning, specifically during the denning period. The findings of this study can be used to reduce the negative effects of existing and planned infrastructure development such as the establishment, operation, or extension of ski resorts, quarries, and other tourism-or transport-related infrastructure, including the planned Tîrgu Mureș-Iași [A8] highway. In the foothills, disturbance by pastoralists, villagers, or domestic dogs constitutes additional disruptions that may be harder to regulate, requiring educational programs in local communities.

Lastly, driven hunts for wild boar (Sus scrofa) should be avoided altogether in prime brown bear denning habitat during the denning period. Driven hunts present the highest chances for bear dens to be encountered. During these hunts, guides methodically scour forest sections for wild boar on foot, driving (chasing) them toward hunters who lie in wait. Such hunts in winter denning areas may result in den abandonment and threat to human safety.

Given the legal status of brown bears as strictly protected in the EU, through the EU Habitats Directive (Council Directive 92/43/EEC), conferring legal protection to areas highly conducive to bear denning will address current and developing threats associated with this key period in the life cycle of brown bears in Romania.

Fig. 1 .

 1 Fig. 1.-Structural components of a typical brown bear den.

Fig. 2 .

 2 Fig. 2.-A map of the study area (foothills in the West, mountains in the East), known brown bear dens (n = 115), and open nests (n = 8) in eastern Transylvania, Romania, 2008-2013 and 2015-2017.

Fig. 3 .

 3 Fig. 3.-Probability of brown bear den presence using data collected from 119 known den sites in eastern Transylvania, Romania. Relative percent contribution of terrain ruggedness to the MaxEnt model is indicated in brackets. Gray area represents the confidence interval for the probability curve among 10 replicates.

  , Eurasian lynx (Lynx lynx-White et al. 2015), and gray wolves (Canis lupus-Ahmadi et al. 2013), and higher terrain ruggedness

Fig. 4 .

 4 Fig. 4.-Potential areas for brown bear denning in eastern Transylvania, Romania, based on MaxEnt algorithm.

Table 2 .

 2 -Measurements for brown bear (Ursus arctos) dens in the mountain (n = 74) and foothill (n = 41) regions of the study area in eastern Transylvania,Romania, 2008[START_REF] Kramer-Schadt | The importance of correcting for sampling bias in MaxEnt species distribution models[END_REF] and 2015-2017. Height (h), width (w), and length (l) values (cm) represent mean ± SD.

		Entrance		Tunnel a			Chamber		Porch b
		h	w	h	w	l	h	w	l	l
	Mountain dens	62.9 ± 23.8	72.9 ± 29.5	72.0 ± 27.7	87.8 ± 32.7	133.7 ± 69.3	82.1 ± 28.6	118.1 ± 36.1	129.7 ± 34.6	171.9 ± 57.9
	Foothill dens	58.3 ± 13.3	82.5 ± 20.3	67.4 ± 15.2	76.5 ± 16.8	112.5 ± 82.4	90.3 ± 29.7	106.5 ± 27.1	129.8 ± 54.1	154.3 ± 62.7

a Only for dens with tunnels (tunnel measurements ≠ 0; 40 dens in the mountains and 11 dens in the foothills). b Only for dens with porches (porch length ≠ 0; 51 dens in the mountains and 32 dens in the foothills).

Table 1 .

 1 -Covariates selected for modeling brown bear den distribution in eastern Transylvania, Romania.

	Covariate ID	Description	Covariate type	References
	TRugg	Terrain ruggedness index (TRI) based on 30-m DEM a	Continuous	Nielsen et al. (2004)
				Whiteman et al. (2017)
	ForestC	Percentage of forest cover within a 495-m rectangle	Continuous	Fernández et al. (2012)
		buffer based on CLC b 2012		Pop et al. (2018)
	EastI	Easterness-aspect index based on 30-m DEM a	Continuous	Craighead and
				Craighead (1972)
				Vroom et al. (1980)
				Schoen et al. (1987)
	NorthI	Northerness-aspect index based on 30-m DEM a	Continuous	Craighead and
				Craighead (1972)
				Vroom et al. (1980)
				Schoen et al. (1987)
	RoadD	Density of paved road within a 500-m circle buffer	Continuous	Swenson et al. (1997)
		based on OSM c (km 2 )		Sahlén et al. (2011)
	ArtifC	Percentage of human infrastructure cover within a	Continuous	Elfström and Swenson (2009)
		495-m rectangle buffer based on CLC b 2012		Fernández et al. (2012)
	OpenC	Percentage of open field cover within a 495-m rec-	Continuous	Elfström et al. (2008)
		tangle buffer based on CLC b 2012

a Digital Elevation Model from the GMES RDA project (EU-DEM). b CORINE Land cover 2012.

c OpenStreetMap project database.

Acknowledgments

We express our gratitude to all who assisted us during fieldwork, most notably game wardens Károly Illyés ("Mureșul" Private Forestry Service), Károly Pál ("Pro Diana" Hunting Association), and Ágoston Pál ("Târnava Mare" Hunting and Angling Association). We thank the foresters of the "Mureșul" Private Forestry Service (Sovata) and to those of the Private Forestry Service from Praid, István Birtalan ("Gheorgheni" Hunting and Angling Association), Béla Lázár and Szilárd Székely ("Cervus" Hunting Association), as well as locals and pastoralists from Neaua. We are grateful to DVM Levente Borka-Vitális (Vets4Wild), Károly Illyés, Károly Pál, and Ágoston Pál for assistance during the capturing and collaring of bears. We thank John Hopkins and two anonymous reviewers for constructive feedback that greatly enhanced the quality of the manuscript. Peter Damerell reviewed and improved earlier versions of the manuscript. The study occurred in the framework of Milvus Group's Brown bear conservation and research program in a model area in Romania. This long-term initiative has received financial support from Bears in Mind (the Netherlands), Bernd Thies Foundation (Switzerland), Columbus Zoo and Aquarium (United States), EuroNatur (Germany), Frankfurt Zoological Society (Germany), International Association for Bear Research and Management (IBA), as well as the Nando Peretti Foundation (Italy). UF was supported by the European Action Scheme for the mobility of University Students (ERASMUS) at the University of Aix Marseille Univ. BC was supported by an Alberta Ingenuity Fund Ph.D. Scholarship, the Alberta Conservation Association, and subsequently a Claude Leon Foundation Postdoctoral Fellowship at the University of Cape Town.

Supplementary Data

Supplementary data are available at Journal of Mammalogy online.

Supplementary Data SD1.-Spatial distribution of the terrain ruggedness index (TRI) covariate selected for modeling brown bear den distribution in eastern Transylvania, Romania.

Supplementary Data SD2.-Spatial distribution of the forest cover (%) covariate selected for modeling brown bear den distribution in eastern Transylvania, Romania.

Supplementary Data SD3.-Spatial distribution of the Northness index covariate selected for modeling brown bear den distribution in eastern Transylvania, Romania. The covariate ranges from -1 (South exposures) to 1 (North exposures).

Supplementary Data SD4.-Spatial distribution of the Eastness index covariate selected for modeling brown bear den distribution in eastern Transylvania, Romania. The covariate ranges from -1 (West exposures) to 1 (East exposures).

Supplementary Data SD5.-Spatial distribution of the paved-road density covariate selected for modeling brown bear den distribution in eastern Transylvania, Romania.

Supplementary Data SD6.-Pairwise Pearson's correlation coefficients (r) between the seven covariates selected for modeling brown bear den distribution in eastern Transylvania, Romania. Covariate details are provided in Table 1.

Supplementary Data SD7.-Structural characteristics of brown bear dens in the mountain versus foothill regions of the study area in eastern Transylvania, Romania, 2008[START_REF] Kramer-Schadt | The importance of correcting for sampling bias in MaxEnt species distribution models[END_REF]and 2015-2017. Statistical comparisons of height (h), width (w), and length (l) are presented for the full data set as well as excluding outliers. Wilcoxon rank sum test statistics (W) as well as probability values (P) are presented. Significant differences (alpha = 0.05) are given in bold.

Supplementary Data SD8.-Ranking of models for brown bear denning habitat distribution in eastern Transylvania, Romania based on ∆AICc. Models were run using MaxEnt and set with and without a bias file as well as linear (L), product (P), and quadratic (Q) features and their combinations.

Supplementary Data SD9.-Model standard deviation for brown bear denning in eastern Transylvania, Romania, based on MaxEnt algorithm, which used data from 119 known den sites.

Supplementary Data SD10.-Probability of brown bear den presence using data collected from 119 known den sites in eastern Transylvania, Romania. Relative percent contributions of four environmental covariates with < 5% individual contributions to the MaxEnt model are indicated in brackets. Gray areas represent the confidence intervals for each probability curve among 10 replicates. Supplementary Data SD11.-Model Multivariate Environmental Similarity Surface (MESS) for brown bear denning in eastern Transylvania, Romania, based on MaxEnt algorithm fitted using 119 occurrences of wintering structures with linear and quadratic features. MESS map shows areas where covariate values are included inside (dark blue) or outside (dark red) the range of covariate values within the area used to fit the model. The map indicates novelty of the environmental conditions within the projected area, and then the degree of extrapolation of the model.