We develop a wide general theory of bilinear bi-parameter singular integrals T . First, we prove a dyadic representation theorem starting from T 1 assumptions and apply it to show many estimates, including L p × L q → L r estimates in the full natural range together with weighted estimates and mixed-norm estimates. Second, we develop commutator decompositions and show estimates in the full range for commutators and iterated commutators, like [b1, T ]1 and [b2, [b1, T ]1]2, where b1 and b2 are little BMO functions. Our proof method can be used to simplify and improve linear commutator proofs, even in the two-weight Bloom setting. We also prove commutator lower bounds by using and developing the recent median method.