We present a general formula describing the joint distribution of two permutation statistics—the peak number and the descent number—over any set of permutations whose quasisymmetric generating function is a symmetric function. Our formula involves a certain kind of plethystic substitution on quasisymmetric generating functions. We apply this result to cyclic permutations, involutions, and derangements, and to give a generating function formula for counting permutations by peaks, descents, and cycle type. Along the way, we recover as special cases results previously derived by Gessel–Reutenauer, Fulman, Fulman–Diaconis–Holmes, and Athanasiadis. This is joint work with Ira Gessel.