We recently carried out a series of DNSs of incompressible turbulence with number of grid points up to 4096(3) on the Earth Simulator (ES). The 4096(3) (ultra) DNS is 64 times the size of the largest DNS before the appearance of the ES. The DNSs consist of two groups; one is with K-max eta = 2 R-lambda = 1 and the Taylor micro scale Reynolds number R-lambda = 167 similar to 1131, while the other is with K-max eta = 2 and R-lambda = 94 similar to 675, where K-max is the maximum wave number and eta the Kolmogorov length scale. A review is given on the one-point and two-point statistics of the velocity derivatives, including the fourth-order rotational invariants I-alpha (alpha = 1...4) proposed by Siggia (1981) and the spectra of energy dissipation epsilon, enstrophy omega(2), and pressure p, obtained by the series of the DNSs on the ES. The series of the DNSs suggest that the spectra of epsilon, omega(2), and p exhibited a wave number range in which the spectra scaled with the wave number k as k(-a). Exponent a for p was in good agreement with the value obtained by assuming the joint probability distribution of the velocity field to be Gaussian, while a values for epsilon and omega(2) were very different from the Gaussian approximation values.