One of the most well-established tools for modeling the brain as a complex system is the functional connectivity network, which examines the correlations between pairs of interacting brain regions. While powerful, the network model is limited by the restriction that only pairwise dependencies are visible and potentially higher-order structures are missed. In this work, we explore how multivariate information theory can reveal higher-order, synergistic dependencies in the human brain. Using the O-information, a measure of whether the structure of a system is redundancy- or synergy-dominated, we show that synergistic subsystems are widespread in the human brain. We provide a mathematical analysis of the O-information to locate it within a larger taxonomy of multivariate complexity measures. We also show the O-information is related to a previously established measure, the Tononi-Sporns-Edelman complexity, and can be understood as an expected difference in integration between system scales. Highly synergistic subsystems typically sit between canonical functional networks, and may serve to integrate those networks. We then use simulated annealing to find maximally synergistic subsystems, finding that such systems typically comprise $\approx$10 brain regions, also recruited from multiple canonical brain systems. Though ubiquitous, highly synergistic subsystems are invisible when considering pairwise functional connectivity, suggesting that higher-order dependencies form a kind of ``shadow structure" that has been unrecognized by established network-based analyses. We assert that higher-order interactions in the brain represent a vast and under-explored space that, made accessible with tools of multivariate information theory, may offer novel scientific insights.