While bending and twist can be treated as independent degrees of freedom for linear DNA molecules, the loop closure constraint introduces a coupling between these variables in circular DNA. We performed Monte Carlo simulations of wormlike rods with both bending and twist rigidity in order to study the coupling between the writhe and twist distributions for various DNA lengths. We find that for sufficiently short DNA, the writhe distribution differs from that of a model with bending energy only. We show that the factorization approximation introduced by previous researchers coincides, within numerical accuracy, with our simulation results, and conclude that the closure constraint is fully accounted for by the White-Fuller relation. Experimental tests of our results for short DNA plasmids are proposed.