In this paper, we study some supercongruences involving multiple harmonic sums by using Bernoulli numbers. Our main theorem generalizes previous results by many different authors and confirms a conjecture by the authors and their collaborators. In the proof, we will need not only the ordinary multiple harmonic sums in which the indices are ordered, but also some variant forms in which the indices can be unordered or partially ordered. It is a crucial fact that the unordered multiple harmonic sums often behave better than the corresponding ordered sums when one considers congruences. We believe these unordered sums will play important roles in other studies in the future.