Here we present a study on three-dimensional shape optimization of flapping wings. A vertical force maximization problem for insect flight is formulated in the domain of unsteady-state, low Reynolds number viscous flow field. Adjoint method, a type of sensitivity analysis based on the variational principle, is modeled by the Navier-Stokes equations and discretized with the finite volume method. With the method, shape optimization of the wings of a model fruit fly in hovering flight is attempted. Rectangles, having both the aspect ratios and the surface areas identical to those of the realistic wings, are taken as the initial shapes. As a result, their sensitivities to the shape modifications are derived.