This study presents an integrated aerodynamic optimization framework for high-altitude long-endurance UAVs operating under low-Reynolds-number conditions, combining surrogate-based and adjoint-based optimization approaches. Surrogate modeling optimized airfoil profiles at wing root, kink, and tip stations, achieving a 3.7% chordwise delay in boundary layer transition onset and 1.3% drag reduction at Ma = 0.5/α = 2°. At Ma = 0.2/α = 2°, lift coefficient increased by 4.6%. Three-dimensional implementation of optimized airfoils enhanced the full-configuration aerodynamic performance, yielding a 1.44 improvement in lift-to-drag ratio at high-speed conditions and 0.05 increase in maximum lift coefficient at low-speed regimes. Adjoint optimization further reduced drag coefficient from 0.0347 to 0.0339 under CL = 0.89 constraint while improving maximum lift-to-drag ratio by 0.85. Comparative flow analysis revealed surrogate-optimized configurations delayed transition onset by 3.4% chord versus baseline, with adjoint-optimized wings exhibiting near-elliptical spanwise lift distribution and 12.7% reduction in crossflow separation areas. Q-criterion vorticity analysis demonstrated adjoint optimization effectively attenuated kinetic energy dissipation zones by 9.3%, particularly in the 60–80% chord trailing edge region, through modified flow control strategies that reduced viscous dissipation effects. These complementary optimization approaches collectively enhanced cruise efficiency through synergistic improvements in flow attachment stability and aerodynamic loss mitigation.