In this paper, we discuss a model order reduction method for multiple-input and multiple-output discrete-time bilinear control systems. Similar to the continuous-time case, we will show that a system can be characterized by a series of generalized transfer functions. This will be achieved by a multivariate Z-transform of kernels corresponding to an explicit solution formula for discrete-time systems. We will further address the problem of generalized tangential interpolation which naturally comes along with this approach. We will introduce a reasonable generalization of the linear H2-norm. Based on this concept, we discuss the choice of interpolation points. Furthermore, an efficient discretization of continuous-time systems is provided. The performance of the proposed method is evaluated by means of some numerical examples and compared with the method of balanced truncation for bilinear systems.