Climate change has profound impacts on regional hydrological characteristics in large unregulated continental river basins (LUCRiBs) such as the Athabaasca River Basin (ARB), Canada. A systematic analysis of these impacts is confronted with many challenges. For instance, the performances of general circulation models (GCMs) vary with many factors, e.g. climate variables, geographic locations, temporal scales, and evaluation measures. Mesoscale atmospheric features can barely be provided by coarse-resolution GCMs. Filling this gap by statistical downscaling is further challenged by redundant computations, resulting from spatial climatic similarities, and the complexities of data uncertainties, nonlinear correspondences, normality prerequisites, and multivariate dependencies. Climatic projection may lack a solid GCM-evaluation foundation and a high spatial resolution. These complexities in downscaling may also exist and be coupled with massive computations in integer optimization in hydrological simulation. Furthermore, an integration of these challenges would decrease the reliability of long-term streamflow forecastings for guiding socio-economic development and eco-environmental conservation over LUCRiBs such as the ARB under climate change. To fill the gap of few effective techniques, an integrated hydro-climatic systems analysis framework is developed and applied to the ARB. This framework includes six modules. (a) The multi-dimensional performances of CMIP5 GCMs and their ensemble are evaluated. (b) The climate over the ARB is classified by recursive dissimilarity and similarity inferences. (c) The spatial resolution of GCM is enhanced by recursive multivariate principal-monotonicity inferential downscaling based on (a) and (b). (d) High-resolution climatic projection under four representative concentration pathways (RCPs) are generated by coupling (a) to (c). (e) The correspondence between climate and streamflow is reproduced by Bayesian principalmonotonicity inference based on (b). (f) Modules (d) and (e) are integrated for streamflow forecasting under climate change. A series of findings are revealed while methodological reliability is verified. For instance, the multi-model ensemble has a relatively high modeling accuracy. The climatic conditions over the ARB are classified into 20 classes based on their dissimilarity and similarity. The overall downscaling accuracies are relatively high for temperature and acceptable for precipitation although varying with multiple factors. At the scale of octo-decades, daily minimum temperature would increase by 1.7, 2.3, 2.1 and 3.0 􀔨, daily maximum temperature by 1.4, 1.8, 1.6 and 2.2 􀔨, and daily total precipitation by 0.03, 0.07, 0.08 and 0.16 mm under RCPs 2.6, 4.5, 6.0 and 8.5, respectively. The approach in module (e) is effective at capturing the temporal variability and the multi-year averages of streamflow and the uncertainties of climatestreamflow correspondences. Streamflow tends to increase at the upper and middle reaches and decline at the lower one. The increments of streamflow would be the highest in March and the decrements would be dominated by less flow in July or Summer. Either RCP scenarios or modeling biases are significant for the temporal variability and trends and are insignificant for the overall magnitudes of streamflow. The methods and findings in this study would be helpful for gaining insights into coupled climatic and hydrological systems over the ARB, evaluating the impacts of climate change, guiding regional socio-economic development and eco-environmental conservation, and promoting developations of more advanced climatic and hydrometeorological systems analysis methods.