Assessment of Heart Rate Variability Thresholds from Incremental Treadmill Tests in Five Cross-Country Skiing Techniques

The assessment of heart rate variability (HRV) thresholds (HRVTs) as an alternative of Ventilatory thresholds (VTs) is a relatively new approach with increasing popularity which has not been conducted in cross-country (XC) skiing yet. The main purpose of the present study was to assess HRVTs in the five main XC skiing-related techniques, double poling (DP), diagonal striding (DS), Nordic walking (NW), V1 skating (V1), and V2 skating (V2).Ten competitive skiers completed these incremental treadmill tests until exhaustion with a minimum of one to two recovery days in between each test. Ventilatory gases, HRV and poling frequencies were measured. The first HRV threshold (HRVT1) was assessed using two time-domain analysis methods, and the second HRV threshold (HRVT2) was assessed using two non-time varying frequency-domain analysis methods. HRVT1 was assessed by plotting the mean successive difference (MSD) and standard deviation (SD) of normalized R-R intervals to workload. HRVT1 was assessed by plotting high frequency power (HFP) and the HFP relative to respiratory sinus arrhythmia (HFPRSA) with workload. HRVTs were named after their methods (HRVT1-SD; HRVT1-MSD; HRVT2-HFP; HRVT2-HFP-RSA). The results showed that the only cases where the proposed HRVTs were good assessors of VTs were the HRVT1-SD of the DS test, the HRVT1-MSD of the DS and V2 tests, and the HRVT2-HFP-RSA of the NW test. The lack of a wider success of the assessment of HRVTs was reasoned to be mostly due to the high entrainment between the breathing and poling frequencies. As secondary finding, a novel Cardiolocomotor coupling mode was observed in the NW test. This new Cardiolocoomtor coupling mode corresponded to the whole bilateral poling cycle instead of corresponding to each poling action as it was reported to the date by the existing literature.


Introduction
Within the last three decades, blood lactate curves and gas exchange measurements from incremental exercise tests have been used for the assessment of endurance capacity and training zones, where two submaximal thresholds have been revealed [1]. However, the assessment of the two submaximal thresholds, being Ventilatory thresholds (VTs) or lactate thresholds (LTs), presents limitations such as laboratory dependence, high costs or an invasive nature [2].
At the beginning of the 21 st century, attempts to overcome these limitations began with the development of a new methodology for the assessment of submaximal thresholds based on Heart Rate Variability (HRV). The assessment of these thresholds, often referred to as HRV thresholds (HRVTs), is based on the links that HRV presents with the autonomous nervous system and respiratory sinus arrhythmia (RSA) [3][4][5][6].
The assessment of the first HRVT (HRVT 1 ) has been conducted with complex time varying frequency analysis methods [7][8][9][10] but also with simple time analysis methods [2,11]. The later ones are based on the virtual disappearance of the N-N (i.e. normalized R-R) interval variations, as seen in the trends of mean successive difference (MSD) and standard deviation (SD) of N-N, which represents vagal withdrawal.
The assessment of the second HRVT (HRVT 2 ) has been conducted with time varying [8][9][10][12][13][14] and non-time varying [13,15] frequency analysis methods. The frequency analysis of HRV decomposes the N-N time-dependent signal into its sinusoidal components, obtaining the so-called power spectral density (PSD). This PSD is then under-divided into low and high frequency (HF) bands. The HF band results from breathing (respiratory sinus arrhythmia or RSA) and mainly vagal cardiovascular control, but it is also related to Cardiolocomotor coupling (LOC) when summit to exercises where the upper body is involved [6,13,15]. The latest studies indicate that the most sensitive methods for HRVT 2 assessment in exercises where the upper body strikes contribute to propulsion are time varying methods that discard the power in the HF band (HFP) relative to the LOC (HFP LOC ) to solely use the HFP relative to the RSA (HFP RSA ) [13,15].
To the best of our knowledge, the assessment of HRVTs has been conducted in varied exercise modes, including exercises using upper body movements, such as front crawl swimming [15] and ski mountaineering [13,14], but never in cross-country (XC) skiing. A characteristic of XC skiing is that, like swimming, it consists of different techniques with different relationships between breathing frequency (BF) and poling frequency (PF) [16,17], which would presumably imply different peak patterns in the HF band. The five main techniques used by XC skiers are Diagonal Striding (DS), Nordic Walking (NW), Double Poling (DP), V1 skating (V1) and V2 skating (V2). All techniques seem to share a tendency for an increase in BF-PF coupling with increasing workload intensities. Moreover, it seems that the PF-BF coupling is the weakest in the DS and NW techniques and strongest in the DP technique [13,[16][17][18][19][20].
The main aim of the present study was to assess HRVTs in the five main techniques used by cross country skiers, and the secondary aim was to evaluate how the Cardiolocomotor interactions affected these assessments. For these purposes, five discontinuous incremental treadmill tests until exhaustion were conducted. Two time analysis methods were selected for the assessment of HRVT 1 . The HRVT 1 obtained from MSD was named HRVT 1-MSD , and the HRVT 1 obtained from SD was named HRVT 1-SD . For the assessment of HRVT 2 two non-time varying frequency analysis methods were used. The HRVT 2 obtained from HFP was named HRVT 2-HFP , and the HRVT 2 obtained from HFP RSA was named HRVT 2-HFP-RSA .

Methods
A group of competitive XC skiers performed five incremental tests until exhaustion in a time frame of two to three weeks with at least one day between tests. The tests were performed in randomized order.

Subjects
Ten competitive national level XC skiers (5 male and 5 female), with ages ranging from 19 to 30 years participated in the study (Table 1). All athletes were healthy and had just completed their XC ski race season. All subjects gave written informed consent agreeing to the terms and conditions of the study, and all possible risks and benefits were explained. The study was approved by the Ethics Committee of the University of Jyväskylä. The amount of subcutaneous fat was estimated from a four site skinfold measurement (m. biceps, m. triceps, os. scapula and, os. crista iliaca) [21].

Procedure
Incremental Tests. NW, DP, DS, V1 skating, and V2 skating were the techniques used for the five incremental tests. All tests were performed on a large motor-driven treadmill (Rodby RL3500E, Rodby Innovations, Vänge, Sweden) located in an indoor laboratory setting. All participants wore a harness that was attached to a rope which hung from a metallic frame in the ceiling above the treadmill for safety. Marwe Classic 800 C and Marwe Skate 610 A roller skis (Hyvinkään Kumi, Hyvinkää, Finland) were used. One Way ski poles (One Way Sport, Vantaa, Finland) were equipped with special rubber tips designed to ensure an optimal grip while roller skiing on the motor-driven treadmill (Biomekanikk AS, Oslo, Norway). Before each incremental test, the subjects performed a 10 to 15 minute warm up with the same technique they were performing in the test on that day. In all protocols, workload increased every three minutes. At the end of every stage, the treadmill was stopped for 10 to 15 seconds for blood lactate samples (not used in this study) from the standing subjects. This small timeframe was included in the three minute stage. The highest 20 second mean values obtained during the tests were defined as peak values. The highest 20 second V O 2 mean value (V O 2 peak) among the five tests was defined as the V O 2 max value for the person if two or more of the following criteria were met: a HR within 10 beats of the age predicted maximal HR (220 -age); a respiratory exchange ratio superior to 1.1; or a plateau in V O 2 with increasing workload [13]. The PF was measured with a chronometer by timing ten poling cycles in the last minute of each stage.
In the DP and the V2 protocols inclination was maintained constant at 3% and 5%, respectively, whereas speed increased 2 km/h with every stage from the initial speeds of 8 km/h for men and 7 km/h for women. In the DS and V1 tests, speed was maintained constant at 10 km/ h for men and 9 km/h for women, whereas inclination increased 2% with every stage from the initial 3% inclination. The NW tests followed a protocol that has been used for more than 30 years in Finland (Table 2) [22].
Assessment of First and Second Ventilatory Thresholds. A portable Oxycon Mobile spiroergometer (Viasys Healthcare GmbH, Hoechberg, Germany) was used for the gas measurements. The spiroergometer's main sampling unit was lightweight and attached to the Assessment of First and Second Heart Rate Variability Thresholds. A Suunto t6 (Suunto Oy, Vantaa, Finland) HR monitor with a sampling frequency of 1000 Hz for the recording of the R-R intervals was used throughout each test. From the HRV recordings, the R-R interval series relative to the last 90 seconds of each step were extracted, as HF oscillations during these periods were considered stationary [15]. The R-R interval series were then automatically filtered with the software Firstbeat Health 3.1.1.0 (Firstbeat Technologies Oy, Jyväskylä, Finland) for the correction of the eventual artefacts and ectopic beats, which has been proven to be fast, accurate and reliable [24]. The normalized R-R (also named as N-N) time series, were then analysed by the software Kubios HRV 2.1 (Biosignal Analysis and Medical Imaging Group BSAMIG, Kuopio, Finland). HRVT 1 was assessed from two time-domain HRV analyses methods [11]. MSD and the SD were graphically plotted against workload. HRVT 1-SD and HRVT 1-MSD were set by visual interpretation at the point where there was no further decrease in the trends of these two parameters [11].
HRVT 2 was assessed with two time-domain HRV analyses methods [15]. The HFP and the HFP RSA were plotted against workload. One and sometimes two increases were noticed in the trends of both variables. HRVT 2-HFP and HRVT 2-HFP-RSA were set by visual interpretation at the point from which the last increase started [15].
The N-N series were first converted into an equidistantly sampled form, and the HR trend and part of the LF oscillations were then removed [15]. The PSD was estimated by an autoregressive model of order 12 or 18, depending on the type of test. An order 12 was applied to the DS and NW tests. However, an autoregressive model order 18 was considered to be more appropriate for the DP, V1, and V2 tests, because unlike the previously studied disciplines [13][14][15], these three disciplines use symmetric bilateral upper body strikes. The LF band was set between 0.04 and 0.15 Hz, and the HF band was set between 0.15 and 2.00 Hz. HFP was computed as the spectral power in the referred HF range, whereas HFP RSA was computed as the power from frequencies ranging from 0.04 Hz to a varying cut-off frequency which corresponded to the borderline between the RSA-related and LOC-related HF components [13,15]. The frequencies of the PSD peaks on HFP RSA and HFP LOC were named as fHF and pfHF, respectively. In the occurrence where HFP RSA and HFP LOC could not be divided, the frequency of the density peak was considered to be both fHF and pfHF. To verify that the HFP RSA and HFP LOC components were indeed related to respiration and locomotion, all of the spectrums (n = 321: 45 maximal tests with 5-10 stages each) were visually inspected and fHF and pfHF were compared to the corresponding BF and PF measured by spiroergometer, and chronometer. The differentiation of the two HF band components in the five different tests occurred as follows. In the HF ranges of spectrums in the V2 tests, two clear PSD peaks emerged, the first peak corresponding to RSA and the second peak corresponding to LOC. In the DS tests, the main peak corresponded to RSA and the last and often non prominent peak's frequency (pfHF) corresponded to the double of PF. In the spectrums of the NW tests, two main peaks emerged. The first peak corresponded to RSA, and the second peak to PF. Lastly, in the vast majority of spectrums obtained from the V1 and DP tests, a single peak emerged in the HF band, preventing to differentiate two HF components. Therefore the decision was made to not split HFP in any of the spectrums of the V1 and DP tests, and the frequencies of their main peaks were decided to be both fHF and pfHF.

Statistical Analysis
A between methods agreement was used to evaluate whether there was an agreement or bias between the variables (i.e. VTs, BF and PF) determined from the reference methods (i.e. timer for PF and Ventilatory gases for the rest of variables) and the corresponding variables assessed from the alternative HRV-related methods (i.e. HRVTs, fHF and pfHF). The analyses included (a) an evaluation of the relationships between parameters using Pearson's r correlation coefficients and linear regression lines, (b) an examination of the level of agreement using 95% limits of agreement according to Bland-Altman, and (c) a comparison of mean values using paired t tests. The agreement analysis between VTs and HRVTs was conducted with HR values because of its practical applicability. All the data are reported as the mean ± SD and the statistical significance was set at p 0.05 for all tests. All statistical analyses were performed using IBM SPSS Statistics 20 software (SPSS Inc, Chicago, IL, USA). The normal distribution of the data was verified by One-Sample Kolmogorov-Smirnov Test and the magnitude of the correlations was assessed according to Hopkins' scale [25].

Agreements between Ventilatory and Heart Rate Variability Thresholds
All VTs except the VT 1 of one subject's DP test were assessed. The R-R interval data from five maximal tests (2 DP tests, 2 V1 tests and 1 V2 test) were excluded due to excess of artefacts, and hence, none of their HRVTs could be assessed. Moreover, HRVT 1-MSD could not be assessed in one of the DS tests because of the abnormal behaviour of its trend. Additionally, the HRVT 2-HFP-RSA was not assessed in the DP and V1 tests because it was not possible to divide the HFP spectrums into two components. Besides the aforementioned cases where HRVTs could not be assessed, the reminders of HRVTs were assessed, including the HRVT 1 corresponding to the VT 1 that could not be assessed by the Ventilatory gas exchange method. A summary of the agreement between VTs and HRVTs can be seen in Table 3.
First Threshold. With regards to the SD method, the NW test was the only test that showed statistical difference between HRVT 1-SD and VT 1 (p = 0.007). These variables where significantly correlated only in the DS (r = 0.77, p = 0.009) and NW (r = 0.68, p = 0.031) tests. The biases show that HRVT 1-SD slightly underestimated VT 1 in the DS and V1 tests (-1 and -2 bpm, respectively), but VT 1 was overestimated in the other tests. The 95% limits of agreement were narrowest in NW and DS (-4/20 and -17/15 bpm, respectively), and was widest in DP (-32/51 bpm). Based on the regression lines, HRVT 1-SD explained 0.1% of the total variability of VT 1 in DP, 59.5% in DS, 46.1% in NW, 8.1% in V1 and 7.0% in V2. Therefore, it seems like the agreement between HRVT 1-SD and VT 1 was best in the DS test ( Fig 1A) since their means were not significantly different, the bias was very small, they were strongly correlated and the regression line explained a good amount of the total variability of VT 1 .
With regards to the MSD method, the DP and V1 tests were the only tests that showed statistical differences between HRVT 1-MSD and VT 1 (p = 0.039 and 0.049, respectively). These variables where significantly correlated in the DP (r = 0.88, p = 0.009), V2 (r = 0.77, p = 0.016) and DS (r = 0.75, p = 0.020) tests. All the biases showed a certain underestimation of HRVT 1-MSD over VT 1 ; the biases where smallest in V2 and DS (-3 and -6 bpm, respectively). The limits of agreement were also narrowest in the V2 and DS tests (-18/13 and -21/10 bpm, respectively). Based on the regression lines, HRVT 1-MSD explained 77.6% of the total variability of VT 1 in DP, 56.5% in DS, 5.8% in NW, 3.6% in V1 and 58.5% in V2. Therefore, it seems like the agreement between HRVT 1-MSD and VT 1 was best in the DS (Fig 1B) and V2 tests, since in both tests the variables were not significantly different, were strongly correlated and the agreement intervals were not that large. The equation of the best fitting linear regression line for the V2 test was y = 1.1007x -19.305.  Second Threshold. With regards to the HFP method, the DS, V1 and V2 tests showed statistical differences between HRVT 2-HFP and VT 2 (p = 0.027, 0.001 and 0.019, respectively). These variables where significantly correlated only in the DS (r = 0.80, p = 0.006) and V2 (r = 0.81, p = 0.008) tests. The smallest bias was found in NW (-4 bpm) and the biggest biases were found in V1 and DS (-18 and 11 bpm, respectively). The 95% limits of agreement were narrowest in V2, NW and V1 (-24/8, -23/14 and -37/1 bpm, respectively), and was widest in DP (-45/28 bpm). Based on the regression lines, HRVT 2-HFP explained 28.6% of the total variability of VT 1 in DP, 63.1% in DS, 19.1% in NW, 7.5% in V1 and 65.5% in V2. Therefore it seems that in none of the tests HRVT 2-HFP was a particularly good assessor of VT 2 . The lack of agreement between HRVT 2-HFP and VT 2 of the NW test is well seen in Fig 2A. With regards to the RSA method, none of the three tests showed statistical differences (p 0.05) between HRVT 2-HFP-RSA and VT 2 . These variables were significantly correlated only in the NW (r = 0.82, p = 0.026) and DS (r = 0.66, p = 0.036) tests. The biases of all three tests were negative, representing an underestimation of HRVT 2-HFP-RSA over VT 2 . The smallest bias was found in NW (-1 bpm) and the biggest in DS (-9 bpm). The limits of agreement was widest in DS (-35/17 bpm) and based on the regression lines, HRVT 2-HFP-RSA explained 44.1% of the total variability of VT 2 in DS, 48.3% in NW and 13.3% in V2. Therefore, it seems like the agreement between HRVT 2-HFP-RSA and VT 2 was close to optimal in the NW test (Fig 2B), since apart from being very strongly correlated and presenting no significant difference, the bias was

fHF-BF and pfHF-PF Agreements
For all tests except for the V2 test, there was no significant difference between the mean of fHF and the mean of the timed BF. These variables were very strongly correlated (r 0.8) in DS, NW and V2, and strongly correlated (r = 0.6-0.8) in DP and V1 (Table 4A). In regards of the agreement between pfHF and the timed PF, the variables were significantly different and uncorrelated (p 0.05) in the DS and V1 tests. Other than this, there were no significant differences in the rest of the tests and the level of correlations were either very strong (NW and V2) or strong (DP) ( Table 4B).

Relationships between Breathing and Poling Frequencies
In the five tests, there was a tendency for an increased coupling between the BF and PF as workload increased and there was always a higher PF than BF during the initial workloads. In the DP test, PF and BF seemed to be quite similar; PF was only significantly higher than BF in the first four stages. However, in the V1 and DS tests the BF trend significantly surpassed the PF trend in a cross-like pattern; significant differences between BF and PF where observed in the initial and last stages. With regards to the NW and V2 tests, substantially higher BF than PF values were presented in all except the last two stages (Fig 3).

Discussion
The present study, the first in assessing HRVTs in XC skiing-related techniques, shows results supporting that the assessment of HRVTs in XC skiing is a difficult task due to the strong involvement of upper body movements. More precisely, the VT 1 was only successfully estimated by the HRVT 1-SD and by the HRVT 1-MSD in the DS test, and by the HRVT 1-MSD in the V2 test. The only good assessor of VT 2 was the HRVT 2-HFP-RSA of the NW test. These assessments support our hypothesis that the VTs would be better replaced by HRVTs in techniques where the involvement of the upper body is moderate (DS and NW) and would by worse replaced when the upper body involvement is highest (DP). Furthermore, the frequency analysis of HRV used in the present study gives an interesting view on how the different poling actions of the studied XC skiing related techniques modulate the HR in different ways, among which a previously undocumented mode of Cardiolocomotor modulation was observed. This new mode of Cardiolocomotor modulation was observed in the power spectrums of the HF bands of the NW test, and was related to half of the rate of arm movements.

Assessment of Ventilatory and Heart Rate Variability Thresholds
Despite the determination of both VT 1 and VT 2 being challenging because the plots did not always clearly show two inflection points, all VTs except the VT 1 of one subject's DP test were determined. These difficulties were expected as a result of a previous V2 study [18] which concluded that the determination of VT 2 from Ventilatory gas analysis must be used cautiously in exercises involving upper body movements. It was explained that the PF entrains the BF, preventing the usual inflection-like behaviour of the Ventilatory variables. With regards to the assessment of HRVTs, taking into account that we hypothesized difficulties for obtaining HRVT 2-HFP-RSA in the DP test, it was not a surprise that HRVT 2-HFP-RSA could not be assessed in DP and V1. This occurred because the strong synchronisation between PF and BF provoked the appearance of single merged peaks in the HF band. However, the reminders of HRVTs were assessed, including the VT 1 which could not be assessed from the ventilatory method.
The HFP RSA and HFP LOC patterns of the studied tests, disregarding the NW test, were not unexpected as they closely corresponded to the frequencies of BF and PF. The appearance of HFP LOC corresponded to the rate of arm movements, which was the same as the PF in V2, V1 and DP, all using the bilateral synchronous polling pattern. In the DS test, HFP LOC also corresponded to the rate of arm movements, but in this particular case it was the same as the double of PF due to DS using a bilateral asynchronous poling (i.e. each poling cycle consists of two poling actions, one for each hand), which was reported in front crawl swimming [15] and ski mountaineering [13]. The same HFP LOC -PF relationship presented in DS was expected to be observed in the NW test, however, this was not the case. The frequency of HFP LOC corresponded to PF, or in other words, to half of the rate of arm movements. As a result the HFP LOC corresponds to the whole bilateral poling cycle instead of each poling action. Thus, during these four techniques using roller skies, the pfHF closely corresponded to the rate of arm movements (i.e 2 PF in DS and PF in V2, V1 and DP), whereas in NW pfHF corresponded to half of the rate of arm movements, "ignoring" the action of half of the poling actions. The correspondence of HFP LOC to half of the rate of arm movements, to the best of our knowledge, is unreported, and it could be considered to be a new Cardiolocomotor coupling mode.
We can speculate that the appearance of this novel LOC mode and the failure to find a peak at the double of PF must have been related, in one way or another, to the distinctive nature that NW is the only discipline that does not use roller skies. On one hand, the absence of a stride phase shortens the stride length, which demands a higher PF for maintaining a given workload.
This implied that, in many spectrums (in 47 out of 87 spectrums) the rate of arm movements was higher than the upper frequency limit of PSD (2 Hz), which would explain the failure to find a peak at the rate of arm movements. On the other hand, the absence of a stride phase most likely also entails lower poling forces due to their propulsion-wise transcendance being reduced. These lower forces could somewhat explain the "missing" peaks in the HF band.

Agreements between Ventilatory and Heart Rate Variability Thresholds
The agreement between VT 1 and the two different equivalents assessed from time-domain analysis were not as good as what were found in a cycle-ergometer [11] and a walking study [2]. The only test where our HRVT 1-SD closely corresponded to VT 1 was during the DS test, and HRVT 1-MSD was only a relatively good assessor of VT 1 in the V2 and DS tests. In the DS test the means of HRVT 1-SD and VT 1 were not significantly different, showed virtually no bias, and were strongly correlated (r = 0.77). However, previous studies [2,11] reported stronger correlations (r = 0.95 and or r = 0.89, respectively). Moreover, the linear regression of our DS test only explained 59.5% of the total variability compared to the 89.6% explanation reported in the walking study [2]. In our V2 and DS tests, the means of HRVT 1-MSD and VT 1 were not significantly different and were strongly correlated (r = 0.77 and 0.75, respectively). However, the cycle ergometer study [11] has again reported a stronger correlation (r = 0.89). Furthermore, the cycle ergometer study showed no bias between the means. On the contrary, the negative biases (-3 and -6 bpm for the V2 and DS tests, respectively) found in our study represented a slight underestimation of HRVT 1-MSD over VT 1 . In the reminder of the tests, both HRVT 1-SD and HRVT 1-MSD failed to successfully correspond to VT 1 .
It is remarkable that, in the cycle ergometer based study [11] where both HRVT 1-SD and HRVT 1-MSD were assessed, the two thresholds corresponded to the same workloads in every case, whereas this was not the case in the present study. Based in our study, it seems that the use of the HRVT 1-MSD is preferable to HRVT 1-SD in XC skiing since the former was a valid assessor of VT 1 in two techniques (V2 and DS) whereas the latter was only a valid assessor in the DS technique.
With regards to HRVT 2 , the agreement between VT 2 and different equivalents assessed from frequency-domain analysis usually did not correspond to the close agreements reported by a swimming study [15] and a ski mountaineering study [13]. Nevertheless, it is important to note that the swimming study [15] compared the HRVT 2 to the second lactate threshold and not to the VT 2 . Our HRVT 2-HFP was not a proper assessor of the VT 2 in either of the tests. However, it seems like HRVT 2-HFP-RSA was a good assessor of VT 2 in the NW test and a fair assessor in the DS test, even though it was not a valid assessor in V2. In the NW test, the means of HRVT 2-HFP-RSA and the VT 2 were not significantly different, they had virtually no bias and they were strongly correlated (r = 0.82). Moreover, the agreement interval was -15 to 14 bpm and the regression line explained 48.3% of the total variability. These results are similar to the results reported in the aforementioned studies [13,15]. In these studies, correlations were also high (r = 0.63 and r = 0.93, respectively), the biases were low (-3 and -1 bpm, respectively), and the agreement intervals were quite narrow (-22 to 14 bpm and -7 to 5 bpm, respectively). In the DS test, the agreement between HRVT 2-HFP-RSA and VT 2 was not as good as in the NW test or in the referred two studies, due to the HRVT 2-HFP-RSA of the DS test considerably underestimating VT 2 (-9 bpm bias). Therefore, HRVT 2-HFP-RSA would probably not be a valid enough alternative of VT 2 in the DS test regarding the applicability in the field, due to its underestimation exceeding the magnitude of the daily variations in HR. A previous study [26] supports this statement, as it concluded a change in submaximal HR of more than 3 bpm could be considered a meaningful change under controlled conditions. Finally, it is also important to note that in our results, similar to the aforementioned literature, when HRVT 2-HFP-RSA was determined, it was a better assessor of VT 2 than HRVT 2-HFP .
A V2 study [18] reported that the assessment of VT 2 in exercises involving upper body movements needs to be used cautiously because there is a high entrainment between the BF and PF. With increasing intensity ventilation shows two inflection points during the regular assessment of VTs, when PF is an added variable, PF alters ventilation and we no longer can observe the inflection points. As a result, the rhythm of PF alters the usual Ventilatory response in which the assessment of VTs (and hence also the assessment of HRVTs) is based on. This could be one reason why the assessment of HRVTs in this study was not always successful. However, some studies [13][14][15] successfully assessed HRVT 2 using the same methods in exercises using upper body locomotion. Finally, as apart from the nature of the exercises, the selected incremental treadmill protocols and the method used for the VTs' assessment were the only differences found in these studies [13,15], it is possible that the selected protocols and the methods used for the VTs' assessment could be improved upon in the present study.

Breathing and Poling Frequencies
With regards to fHF-BF and pfHF-PF agreements, BF and PF assessed by HRV (i.e fHF and pfHF) and their respective equivalents determined conventionally did not present any high variances. If any, the HRV measurements tended to slightly overestimate the reference measurements. This was expected from the literature [13,15,27], and thus, the study was successful in regards to this matter.
With regards to the relationships between breathing and poling frequencies, all five tests had a tendency to show an increased BF-PF coupling with increased workload and a higher PF than BF during the initial workloads. This is in line with the previous literature [18]. With regards to the five specific tests, three kinds of BF-PF relationships were observed: 1) In the DP test, BF and PF were synchronised (1:1 BF to PF ratio) almost from the very beginning; 2) the V1 and DS tests showed statistical differences in early and late stages as the BF trend surpassed the PF trend in a cross-like pattern; and 3) the NW and V2 tests presented substantially different BF and PF values except in the last two stages. All results agree with the literature [13,[16][17][18][19][20] except for the significantly different mean values observed in the final stages of the V1 and DS tests that lead to the cross-like pattern. This cross-like pattern phenomenon, to the best of our knowledge, was unreported to date. Moreover, our initial hypothesis stating that there would be an especially high BF to PF ratio for the DS test and a lower BF-PF entrainment for the DS and NW tests was confirmed. Nevertheless, apart from the DP test, the V1 test also showed a 1:1 BF to PF ratio predominance among the subjects, and in the V2 test, a high BF-PF correlation was also confirmed. In the V2 test, most of the subjects adopted 1:1 and/or 1:2 BF to PF ratios as reported in the previous literature [16,18].

Conclusion
The assessment of HRVTs is a field of research worthwhile to study as it aims to offer an economical, non-laboratory dependant and non-invasive alternative to the conventional Ventilatory and lactate thresholds. Nevertheless, this paper shows that the assessment of HRVTs, when applied to XC skiing-related techniques, is a difficult task that requires expertise due to the poling actions altering the usual Ventilatory response in which the assessment of HRVTs is based off of. Moreover, this paper also shows that the assessment of HRVTs in XC skiing requires the selection of an appropriate method, which might not always ensure the successful determination of HRVTs. The only cases where this study successfully assessed VT 1 from HRV were the two time-domain HRV methods in the DS test and the MSD-based method in the V2 test, whereas the VT 2 was only succesfully assessed in the NW test with the HFP RSA -based frequency-domain HRV method. Thus, one could encourage the use of the time-domain HRV methods in the aforementioned cases where they were proven successful, because they are cheap and fairly simple, although the applicability of the proposed frequency-domain analysis methods for practical settings are questioned due to being more complex, and hence, knowledge-and time-demanding.
We would propose to single out specific XC skiing protocols, using the DS technique to determine HRVT 1 and the NW test for the assessment of HRVT 2 , rather than using all five techniques. Moreover, its applicability is yet to be evaluated in training routines. In any case, further research is needed to prove the validity and reliability of the proposed methods for the assessment of HRVTs in training routines, especially because there seem to be quite wide agreement interval for individual purposes. Further research is also needed for finding alternative HRV-based methods that will successfully assess the VTs in the cases where the present study failed. We suggest the use of nonstop incremental treadmill test protocols, with breathby-breath analysis methods and time varying frequency-domain analysis of HRV for the assessment of the VT 2 . This would be even more time-demanding than the non-time varying frequency domain method used in this study, but it has already proven to be more efficient in ski-mountaineering [13,14].