Activator-inhibitor FitzHugh-Nagumo (FHN) equation is an example for reaction-diffusion equations with skew-gradient structure. We discretize the FHN equation using symmetric interior penalty discontinuous Galerkin (SIPG) method in space and average vector field (AVF) method in time. The AVF method is a geometric integrator, i.e. it preserves the energy of the Hamiltonian systems and energy dissipation of the gradient systems. In this work, we show that the fully discrete energy of the FHN equation satisfies the mini-maximizer property of the continuous energy for the skew-gradient systems. We present numerical results with traveling fronts and pulses for one dimensional, two coupled FHN equations and three coupled FHN equations with one activator and two inhibitors in skew-gradient form. Turing patterns are computed for fully discretized two dimensional FHN equation in the form of spots and labyrinths. Because the computation of the Turing patterns is time consuming for different parameters, we applied model order reduction with the proper orthogonal decomposition (POD). The nonlinear term in the reduced equations is computed using the discrete empirical interpolation (DEIM) with SIPG discretization. Due to the local nature of the discontinuous Galerkin method, the nonlinear terms can be computed more efficiently than for the continuous finite elements. The reduced solutions are very close to the fully discretized ones. The efficiency and accuracy of the POD and POD-DEIM reduced solutions are shown for the labyrinth-like patterns.