The relationships between open-shell characters (y) and second hyperpolarizabilities (γ) of hexagonal graphene nanoflakes including boron nitride domains (BN─HGNFs) are investigated using the long-range corrected spin-unrestricted density functional theory method. The y values are found to alternately decrease and increase with the BN-domain size in a similar manner to the antidot-size dependence of y in antidot HGNFs. It turns out that the variation of y significantly affects their γ values in accordance with the y─γ correlation: the largest γ of the intermediate open-shell BN─HGNFs is about 3.5 times as large as that of the closed-shell analog.