This work presents Bayesian network-informed conditional random forests (BNICRF), a novel multiresponse classification technique able to downscale the joint probability distribution of precipitation occurrence at a set of geographical locations of interest from large-scale predictors. BNICRFs combine a Bayesian network, which provides a description of the multisite/spatial dependence structure of the target locations; with an ensemble of conditional random forests, which extract the relevant information from the large-scale predictors to produce accurate probabilistic predictions. This work extends two previous articles that explored Bayesian networks (Legasa & Gutiérrez, 2020) and random forests (Legasa et al., 2022) in the context of statistical climate downscaling. Building on the experimental and validation framework proposed in the experiment 1c of the COST action VALUE (the largest, most exhaustive intercomparison study of statistical downscaling methods to date) and under the perfect prognosis approach, we thoroughly assess the performance of the proposed methodology focusing on its capability to reproduce several multisite and single-site statistics of the observed series. BNICRFs accurately capture the relevant spatial relationships while keeping the same single-site predictive performance than both single-site random forests and a Generalized Linear Model that exhibited a good performance in VALUE. Moreover, we compare BNICRFs against a robust multisite methodology proposed in Chandler (2020), obtaining better predictive capability while keeping similar spatial performance. Since the assessed models still underestimate autocorrelation, we also propose a straightforward extension of BNICRFs that incorporates the temporal structure at each location of interest able to produce both temporarily and spatially realistic precipitation occurrence fields